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A COMPUTER PROGRAM FOR AUTOMATED FLUTTER 
SOLUTION AND MATCHED-POINT DETERMINATION 

By Kumar G. Bhatia* 

Langley Research Center 

SUMMARY 

The use of a digital computer program (MATCH) for automated determination of 
the flutter velocity and the matched -point flutter density is described. The program is 
based on the use of the modified Laguerre iteration formula to converge to a flutter 
crossing or a matched-point density. 

A general description of the computer program is included and the purpose of all 
subroutines used is stated. The input required by the program and various input options 
are detailed, and the output description is presented. The program can solve flutter 
equations formulated with up to 12 vibration modes and obtain flutter solutions for up to 
10 air densities. The program usage is illustrated by a sample run, and the FORTRAN 
program listing is included. 


INTRODUCTION 

An automated method for determining the flutter velocity and the matched-point 
flutter density is described in reference 1 which contains the theoretical development of 
the method and outlines the computational steps necessary to implement the method on 
a digital computer. However, reference 1 does not contain detailed information about 
the computer program MATCH that was developed to implement the method. The purpose 
of this report is to serve as a user’s manual for this computer program. The basic equa- 
tions used in the computer program are repeated from reference 1, and the general 
program organization is described. The purpose of all the subroutines used is stated, 
and flow diagrams for the two main subprograms are included. The program input and 
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output are described, and a sample run of the program is included in appendix A. The 
FORTRAN program listing and the Langley library subprograms used by MATCH are 
described in appendixes B and C, respectively. 

The present report relies on reference 1, but this report contains complete infor- 
mation regarding the use of the computer program. It is, however, recommended that 
reference 1 be used in conjunction with this report for a complete understanding of the 
theoretical basis of the procedure implemented, 

SYMBOLS 

[A] nondime ns ional aerodynamic matrix (see eq. (4)) 

[AI] = 4jt(BR) 3 (||) 2 (l) 2 [A] 

[AF] = p[AI] 

A s airspeed 

BR reference chord length 

F " Vf - A s 

\G\ vector of damping functions (see eq. (5)) 

lGl{ first partial derivative of {g} with respect to i 

A 

■(G2} second partial derivative of {g} with respect to ^ 

K 

[i] identity matrix 

\ 

IOK current value of i (see eq. (6)) 

k reduced frequency 
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NM 

number of modes 


»— i 

vector of predicted values of ^ corresponding to flutter crossings 


[SK] 

symmetric structural stiffness matrix 


[SM] 

symmetric structural inertia matrix 


SS 

semi span 


jUmi 

eigenvector (see eq. (1)) 


V 

velocity 


Vf 

lowest flutter velocity for an air density 


j V m[ 

associated eigenvector (see eq. (2)) 


v f 

flutter velocity 


, x >y 

Cartesian coordinates 



eigenvalue (see eqs. (1) and (2)) 


i q ‘ 
Id? 

n 



p 

air density 


<>o 

air density at sea level 


CO 

harmonic frequency 


Superscripts: 


R,I 

denote real and imaginary parts of a complex number, respectively 
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T denotes a matrix transpose 

Subscript: 

m denotes the mode number 

Subscripts following a parenthesis denote derivatives. 

GENERAL DESCRIPTION OF THE COMPUTER PROGRAM 

The basic equations used to implement the procedure for the flutter solution and 
determination of the matched-point flutter condition and the general organization of the 
computer program MATCH are described in this section. A matched-point flutter condi- 
tion is obtained when the flutter velocity, air density, and Mach number are consistent for 
standard atmospheric conditions. The aerodynamic matrices are generated external to 
the present program and are required as input to MATCH. These matrices are calculated 
for a given structural configuration and a fixed Mach number. The reduced frequency 
range of interest is selected, and the aerodynamic matrices are evaluated at discrete 
values of the reduced frequency within the selected range. The Mach number is held 
fixed in the program, and therefore the same set of aerodynamic matrices is used. 

The program can be used to obtain a flutter solution at one or more air densities 
or to determine a matched-point density. For a flutter solution at a specified air density, 
an initial value of the inverse of the reduced frequency is input to start the iteration pro- 
cedure, and the program will automatically determine the velocities at which the damping 
becomes zero, if any, within the range of reduced frequency for which the aerodynamic 
matrices have been input. If a matched-point density is desired, an initial air density 
and an inverse of the reduced frequency are input into the program. The program deter- 
mines the lowest flutter velocity for the input density. This flutter velocity will, in 
general, not be the same as the airspeed corresponding to the input density and the fixed 
Mach number. A new air density is predicted to yield the matched-point flutter condition, 
and the lowest flutter velocity for the predicted air density is determined for comparison 
with the airspeed (at the predicted density). This procedure is repeated until an air 
density is determined where the lowest flutter velocity is within a specified tolerance of 
the airspeed. 
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The program is dimensioned for a maximum of 12 modes and 10 air densities, 
that is, the structural and aerodynamic matrices can be up to (12 X 12), and flutter solu- 
tions for up to 10 air densities may be obtained during one run. The program does not 
provide a rigid-body mode capability, but it is possible to extend the program to include 
rigid -body modes. The program requires a field length of about 46 000 octal storage 
locations plus the field length required by the loader. 

Equations Required To Implement the Procedure 

The basic equations to implement the flutter solution procedure and to determine 
the matched-point flutter density are stated in their final form. The derivation of these 
equations is given in reference 1 and is not repeated here. 

The characteristic flutter equation is expressed as an eigenvalue problem in 
matrix form by 

[SK]-1 [[SM] + [AF]] - |u m | = M (m = 1, . . NM) (1) 


where M m and |u m [ are the complex eigenvalues and eigenvectors, respectively. The 
associated eigenvectors jv m ( are determined from the following equation: 


[SK]-1 [[SM] 


+ [AF]T] - Mm [I[ 



jo| 


(m = 1, . . NM) (2) 


where 


[AF] = p[AI] (3a) 

[AI| = 4 7r (BR)3 (||) 2 (e) 2 [A] (3b) 

and the elements of [A] are nondime ns ional. Each element Ay of matrix [A] is defined by 






S 


APj(x, y) dx 


dy 

(SS/BR) 


(4) 


5 



where h^x, y) is the displacement in the ith vibration mode, andApj(x, y) is the aero- 
dynamic pressure over the lifting surface S induced by the downwash associated with 
simple harmonic motion in the jth vibration mode. 


A flutter solution is obtained (for an assumed density p) when the imaginary part 
of one of the eigenvalues of equation (1) (or eq. (2)) is zero. A damping function G(M) 
is defined for each eignevalue p m and is given by 


G(M) 


Mm 



(m = 1, . . ., NM; M - 1, . . NM) (5) 


where 


Mm 


= M m 
h m 



Thus, a flutter solution is obtained when one of the damping functions is zero and the 
corresponding frequency is real (m^i >0 )- Each G (M) is regarded as a function of the 
inverse of reduced frequency 

A modified Laguerre iteration scheme is used to predict a value of ^ for which 
the damping function would be zero and the slope of the curve for damping as a function 
of is positive. The modified Laguerre formula used to predict a zero of G(M) is 


HFI(M) = IOK - - - (6) 

|/[G1(M)] 2 - [G(M)] [G2(M>] 

where IOK is the current value of ^ and RFI(M) is the predicted value of ~ corresponding 
to G(M) = 0. The first derivative (G1(M)) and the second derivative (G2(M)) of G(M) with 
respect to ^ are evaluated from the following expressions: 

("k)l- G < M >f£)l 

G1(M) ^ — (7) 

m 

and 


G2(M) = 


fm)l 1 ' G(M) ("S)l 1 - 

kk kk 



2G1 <M)«)l 

k 


( 8 ) 
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The expressions for^ m j 1 and (^ m ) l i are given by equations (9) and (10) of 
reference 1. ^ ^ ^ 

The flutter solution is determined by using equations (1) to (8), and there may be 
several flutter crossings where one of the damping functions is zero. The flutter solution 
consists of the values of and J 2 and v f for the flutter mode (eigenvalue for which the 
damping is zero) at each crossing. The lowest flutter velocity Vj thus obtained will, in 
general, not be consistent with the airspeed A s for the assumed Mach number (fixed) and 
the assumed air density determined from the standard atmosphere (ref. 2). An iterative 
scheme similar to that used for the flutter solution is used to predict an air density at 
which the lowest flutter velocity and the relevant airspeed will be nearly the same. A 
function F is defined as 


F = F(£ ) = V f - A s 


( 9 ) 


where £ = > and Vf and A s are also regarded as functions of £ . It is apparent that 

a zero of F will yield a matched-point density. The predicted zero of F, that is, £ p , is 
determined from the Laguerre formula 


2F 


£o = £ 


(F) e + sgn[(F) ? ]{[(F) £ ] 2 - 2F(F)jJ 


1/2 


( 10 ) 


where (F)^ and (F)^, respectively, are evaluated by using equations (18a) to (20f) of 
reference 1. A flutter solution is again obtained for the air density corresponding to £ p , 
and the value of F is determined. If F is within some acceptable tolerance, then the 
iteration is terminated; if F is not within the tolerance, the whole procedure is repeated. 


Organization of Program MATCH 

The program MATCH is divided into the two major subprograms LEFCROS and 
CROSMAT. LEFCROS is the subprogram which controls the basic flutter solution capa- 
bility, and CROSMAT controls the determination of the matched point. Both of these 
subprograms call various other subprograms, and since flutter solutions are required 
as a part of the matched -point search CROSMAT calls LEFCROS. Simplified flow dia- 
grams of subprograms CROSMAT and LEFCROS are presented and the various sub- 
routines called by these two subprograms are described subsequently in this section. 
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Simplified Flow Diagram of Subprogram CHOSMAT 


IMATCH = 1 




CALL LEFCROS 


CALL SOUND 


F = V f - A s 


PER = F x 


100/V f 


PER % PERF? 


Legend: 



PERF 


indicates call to a subroutine 


convergence tolerance for 
matched-point flutter 
velocity 


CALL DERVDEN 


H = [(F)!] 2 - 2F(F)|! 


,H g 0? 


Predict £p 


p - «p 2 Ao 


Predict 1A for 
lowest flutter 
velocity crossing 


CALL FOMATCH 


IMATCH = IMATCH + 1 















Simplified Flow Diagram of Subprogram LEFCROS 


Legend; 


Entry FOMATCH 


indicates call tu a subprogram 

number of densities for which flutter 
solution is required 
number of flutter crossings to be 
searched for each density 
maximum number of iterations 
per crossing 

guess for 1/k at a crossing 
minimum and maximum values oi l/k, 
respectively, for which aerodynamic 
information is available 
a matched-point flutter solution is 
not required 
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The aerodynamic matrices and their derivatives are stored in the program on a 
random- access file for easy retrieval during program execution. This aerodynamic 
information is furnished as input to the program by the user. The random-access file 
is generated in subprogram HAND AX which is called from subprogram LEFCROS. The 
required aerodynamic matrix and its first two derivatives are retrieved from the random- 
access file by calls from LEFCROS to subprogram GETAERO and entry point GETDAER. 
These and other subprograms called from CROSMAT and LEFCROS are briefly described. 


Subprogram 


Description 


SOUND Determines the airspeed and its first two derivatives with respect to $ 
between geometric altitudes of -5000 meters (£ = 0.7964651669) to 
20 000 meters {£ = 3.711884976). The airspeed is expressed as a second- 
order polynomial in £ for altitudes between -5000 meters and 11 100 meters, 
and as a constant between 11 100 meters and 20 000 meters. This functional 
representation is based on data from U.S. Standard Atmosphere, 1962 (ref. 2). 

DERVDEN Evaluates the first two derivatives of the flutter velocity with respect to £, 

It also determines the first two derivatives of the reduced frequency and the 
flutter frequency squared with respect to £. It calls subprogram T MM PROD 
to evaluate a matrix triple product. 


FOMATCH 

MATINV 


This is an entry point in LEFCROS, and is called from CROSMAT. 

Langley library subroutine used for determining the inverse of the stiffness 
matrix when the matrix is not diagonal. (See appendix B.) 


RAND AX 


Called only once to read nondimens ional aerodynamic matrices from a disk 
file, convert them to appropriate dimensional form, and write them on a ran- 
dom access file. It uses computer-system-dependent Control Data subrou- 
tines OPENMS and WRITMS at Langley Research Center. (See appendix B.) 


GETAERO Called to retrieve the aerodynamic matrices corresponding to a value of the 
inverse of reduced frequency from the random access file generated by sub- 
program RAND AX. It uses computer -system -dependent Control Data 
subroutine READ MS. (See appendix B. ) 


EIGSOL 


Called from LEFCROS to determine the eigenvalues, eigenvectors, and 
associated eigenvectors by solving equations (1) and (2). It calls Langley 


10 



Description 

library subroutine EECM to solve these equations (see appendix B), and 
subprogram T MM PROD to evaluate triple matrix products. 

This is an entry point in GETAERO and is called from LEFCROS to retrieve 
the derivatives of the aerodynamic matrix from the random access file. 

Called from LEFCROS to evaluate the first two derivatives of the inverse of 
flutter frequency squared with respect to the inverse of reduced frequency. 

It calls subprogram T MM PROD to evaluate the matrix triple products 
required. 

Called from LEFCROS to calculate the predicted values of the inverse of 
reduced frequency corresponding to zero damping crossings by using 
equation (6), and arranging the predicted values in ascending order. It 
calls subprogram DAMPAR to calculate damping from each eigenvalue by 
equation (5). 

INPUT AND OUTPUT DESCRIPTION 
Input 

The input required by the computer program is described. There are two types 
of input to the program: 

(1) Aerodynamic matrices through a disk file (tape 4) 

(2) Namelist input 

Description of tape 4. - All the aerodynamic matrices are in nondimens ional form. 
These matrices must be generated by the user and provided as input to the program on a 
disk file (tape 4) in a format and arrangement that is compatible with the program read 
operations described. Tape 4 is rewound in the program and all information is read in 
binary. 

The first read statement executed is 
READ(4) NK, MACH, NM 

where NK is the number of reduced frequencies for which aerodynamic matrices are on 
tape 4 (NK ^ 1600), MACH is the Mach number at which the aerodynamic matrices have 


Subprogram 

GETDAER 

DERF 

LEGROOT 
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been calculated, and NM is the number of modes defining the size of the aerodynamic 
matrices. The next 3NK read operations are described by the following three read state- 
ments executed NK times: 

HEAD (4) RF,X,((A r (L,M),L=1,NM),M=1,NM),((A i (L,M),L=1 ? NM),M=1,NM) 

READ (4) ((DA r (L, M), L=l, NM), M=l, NM), ((DA X (L, M), L=1,NM), M=l, NM) 

READ(4) ((SDA R (L, M), L=l, NM), M=l, NM), ((SDA*(L, M), L=l, NM), M=l, NM) 

where 

NM number of modes, ^12 

RF reduced frequency for which the six aerodynamic matrices have been 

calculated 

X a dummy scalar (real), not used in the program 

A r real part of (NM X NM) aerodynamic matrix defined by equation (4) 

Al imaginary part of (NM X NM) aerodynamic matrix defined by equation (4) 

DA R first partial derivative of A R with respect to reduced frequency 

DA* first partial derivative of A* with respect to reduced frequency 

SDA R second partial derivative of A R with respect to reduced frequency 

SDA* second partial derivative of A* with respect to reduced frequency 

It is required that the aerodynamic matrices be on tape 4 for increasing values 

of the inverse of reduced frequency ~ and at a constant increment of ^ . For example, 

if the first value of ~ for which the aerodynamic matrices are on tape 4 is RFIL 

(RF^ = l/RFIL), the second value of ^ must be RFI 2 = RFIL + DEL (RF 2 = l/RF^), and 

the last value of i must be RFIR = RFIL + (NK-1) DEL (RF Nk = l/RFIR) where DEL is 

the constant increment in ^ . 

k 

Description of namelist input. - The following two namelists are read from the input 
file in the order presented, 

(1) NAMELIST/NAMATCH/PERF, MAXMAT, MACH, ITROPO, IMATCH 
REFSLD, UNITL 

(2) NAMELIST/NAM1/SK, SM, L STIFF, SS, BR, NM, RFIL, RFIR, DEL, NROOT, 
NITMAX, ND, RHO, R FI MIN, IPRT, IOPT 
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The dimensional parameters in the namelist statements determine the force (for example, 
newtons, pounds, etc. ) and length (meters, feet, etc. ) units with which the program 
operates; the unit of time used is seconds. The user must therefore prepare the name- 
list input to be consistent with any desired force and length units. The definitions of the 
various namelist input parameters in NA MATCH are 


PERF 


Nondimens ional convergence tolerance for matched-point flutter solution. 


The program will terminate when 
Not required if IMATCH = 0. 


1 - 


Airspeed 


Lowest flutter velocity 


x 100 s i erf. 


MAXMAT Maximum number of iterations permitted for the matched -point density 
search. Not required if IMATCH = 0. 

MACH Mach number (real variable) for which the aerodynamic matrices have been 
calculated. Not required if IMATCH = 0. 


ITROPO Defines the initial air density for the matched-point search if 
= 0, initial density = sea-level density 
= 1, initial density = density at altitude of 11 100 meters 
= -1, initial density = RHO(l) from input for namelist NAM1. 
Not required if IMATCH = 0. 


IMATCH 

REFSLD 

UNITL 


If IMATCH = 0, flutter solutions for densities in namelist NAMl are required. 
If IMATCH & 0, a matched -point flutter solution is required. 

Reference sea -level density in ffQ r £fL-s . ep 2 un it s . Not required if IMATCH = 0. 

(Length) 4 

Ratio of the number of length units selected to 1 foot. Not required if the 
length units selected are feet. 


The definitions of the various namelist input parameters in NAMl are 


SK 


SM 

L STIFF 


Force 


Symmetric structural stiffness matrix ^ Length unitsj or symmetric structural 
flexibility matrix unitsj, (NM X NM). 

Symmetric structural inertia matrix ^ ^ength ' c2 ) > ( NM x NM). 


If = 0, SK is diagonal stiffness matrix. 

If = +1, SK is nondiagonal symmetric stiffness matrix. 

If = -1, SK is flexibility matrix, and may or may not be diagonal. 
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ss 

BR 

NM 

RFIL 

RFIR 

DEL 

NROOT 

NITMAX 

ND 

RHO 

RFIMIN 

IPRT 


Semispan (or reference length) used to generate the aerodynamic matrices 
(Length units). 

Reference chord used to generate the aerodynamic matrices (Length units). 

Number of vibration modes used to generate aerodynamic matrices. 
Maximum value of NM is 12. 

Inverse of reduced frequency ( nondimens ional) corresponding to first value 
of reduced frequency for which aerodynamic matrices are on tape 4. 

Inverse of reduced frequency (nondime ns ional) corresponding to last value 
of reduced frequency for which aerodynamic matrices are on tape 4. 

Constant increment of the inverse of reduced frequency (nondime ns ional) 
at which the aerodynamic matrices are on tape 4, for example, the first 
set of matrices are for ^ = RFIL, the second set for ^ = RFIL + DEL, etc. 

Number of flutter crossings desired. If the program cannot determine all 
the NROOT crossings within the selected range of RFIL to RFIR, it will 
continue with execution of the next task, if any. 


Maximum number of iterations per crossing allowed for convergence. If a 
crossing cannot be determined in NITMAX iteration, the execution will be 
terminated. Suggested value 5. 


Number of densities for which a flutter solution is required (is ND s 10). 
If IMATCH* 0 in namelist NAMATCH, then ND should be input as 1. 


One -dimensional array of input densities 


/Force -(Second)^ 
\ (Length) ^ 


unitsj 


for which 


flutter solutions are required. If IMATCH ^0 and ITROPO = -1 in namelist 
NAMATCH, then RHO(l) is the initial density for the matched-point density 
search. If IMATCH^ 0 and ITROPO =£ -1, no input is required for RHO. 


Initial guess for inverse of reduced frequency to start search for first flutter 
crossing. Experience with the program indicates that the convergence from 
a value of RFIMIN which is higher than the inverse of reduced frequency for 
the (actual) first crossing is faster than that from a RFIMIN which is lower. 


Determines amount of output printed by program. It is nominally set to zero 
within the program, and if nominal output is required, then it can be omitted 


14 



from the namelist input. This procedure will be discussed further when 
the program output is described. 

10 PT Unused parameter, not required. 

Output 

The program output is described in this section. The program output consists of 
two categories: 

(1) An output summary on a coded (BCD) disk file (tape 7) which may be routed 
for printing. 

(2) Output file containing either a nominal printout (IPRT = 0 specified by input) 
or a detailed printout (IPRT = 1 or 2). 

The first category of the output is described first and is followed by the second category. 
The output is in all cases in units consistent with those used for the program input. 

The output summary on tape 7 includes the following: 

(1) Air density and the initial value of i for each air density at which a flutter 
solution is determined. 

(2) Root number, flutter velocity, the inverse of reduced frequency, and the total 
number of iterations required for each flutter crossing determined. 

(3) Iteration number, air density, square root of sea -level density/air density, 

lowest flutter velocity, airspeed, and ^ - l qw6 st ^fi^tt e ^ ve lo c ityj x 1^0 * or eac ^ 
matched -point iteration, if matched -point density search is executed. 

(4) Informative messages: 

(a) "FOUND NR ROOTS, RFI FOR THE NEXT ROOT PREDICTED = X, 

IS BEYOND RANGE. " This message is printed when the predicted 
inverse of reduced frequency (X) for the NRth crossing (NR s NROOT) 
is not within the range RFIL to RFIR. 

(b) "RFI PREDICTED FOR THE NEXT ROOT “ . . DIFFERENCE FROM RFI 
FOR PREVIOUS ROOT LESS THAN DEL/2. 0." This message is printed 

to inform the user that the next flutter crossing is within DEL/2. 0 of 
^ for the previous flutter crossing; therefore, the next crossing is taken 
to be at the same value of ^ as the previous flutter crossing. 
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(5) Various messages explaining abnormal termination: 


(a) "MATCH-POINT ITERATION DID NOT CONVERGE IN MAXMAT 
ITERATIONS. M 

(b) "ARGUMENT OF RADICAL IN LAGUERRE = X, ITERATION NO. *. . 
DENSITY * . . VEL = . . SPEED OF SOUND * MACH = . . ." This mes- 
sage is printed out when X is negative during a matched -point density search. 

(c) "ARGUMENT OF RADICAL IN LAGUERRE ITERATE IS NEGATIVE FOR 
NM MODES. " This message is printed out when a real value for predicted 
inverse of reduced frequency for a flutter crossing (from eq. (6)) cannot 
be obtained for any of the NM modes. 

(d) "RFI * . . IS OUTSIDE THE RANGE OF VALUES." This message is 
printed when the initial value of the inverse of reduced frequency input 
in the program is outside the range RFIL to RFIR. 

(e) "PROGRAM TERMINATED, COULD NOT FIND ROOT NO. NR IN NITMAX 
ITERATIONS. " 

(f) "NUMBER OF EIGENVALUES COMPUTED M. " This message is printed 
out from subprogram EIGSOL when during eigensolution, convergence is 
obtained for only M < NM eigenvalues. 

The second category of the output depends on the value for IPRT (0, 1 or 2). In 
all cases, the printout described for tape 7 is included. The output for IPRT = 0 is 
described by stating the additional output relative to printout on tape 7, output for IPRT = 1 
is described by stating the additional output relative to IPRT = 0, and the output for 
IPRT = 2 is similarly described. 

For IPRT = 0, the following information is printed in addition to the information 
written on tape 7: 

(1) Printout of the two namelists. 

(2) Eigenvalues and the predicted values of the inverse of reduced frequency (RFI) 

in increasing order of magnitude, for flutter crossings at each iteration. RFI = 1000, 0000 
indicates that a real value for the inverse of reduced frequency corresponding to a flutter 
crossing could not be predicted for that mode. RFI = 3000.0000 indicates that the real 
part of the eigenvalue corresponding to that mode was negative. 
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(3) Flutter eigenvalue number, eigenvalues, correspondence of the predicted 
crossings (which are arranged in increasing order of magnitude) to eigenvalues (which 
are obtained (and printed) in the decreasing order of their absolute values from the eigen - 
solution), root number, flutter velocity, the inverse of reduced frequency at the crossing, 
and the total number of iterations required for convergence for each flutter crossing 
determined. The number of iterations for convergence include the last iteration for the 
convergence check. 

(4) The inverse of reduced frequency, the first and second derivatives of reduced 
frequency, flutter frequency squared, and flutter velocity with respect to £ , for each 
match-point iteration. 

(5) Predicted values for £ and the inverse of reduced frequency for a matched-point 
solution for each match-point iteration. 

For IPRT = 1, damping, the first and second derivatives of damping with respect 
to ^ , argument of the square root in equation (6), and the predicted crossing are printed 
for each mode during every iteration for a flutter solution. If IPRT = 2, the eigenvectors 
and associated eigenvectors, and the eigenvalue derivatives are printed during every 
iteration for a flutter solution. 


CONCLUDING REMARKS 

A digital computer program MATCH for automated determination of the flutter 
velocity and the matched-point flutter density has been described. The program was 
based on the use of the modified Laguerre iteration formula to converge to a flutter 
crossing or a matched-point density. 

A general description of the computer program and the related subroutines has been 
included. Detailed descriptions of the output, input, and input options have been presented. 
The program can solve flutter equations formulated with up to 12 vibration modes and can 
obtain flutter solutions for up to 10 air densities. Use of the program is illustrated with 
a sample run and the FORTRAN listing is included. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., September 10, 1973. 
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APPENDIX A 


SAMPLE RUN OF PROGRAM MATCH 

The input and output for a sample program run are presented in this appendix in 
order to illustrate the application of the program. The units used in this sample run are 
pounds and inches; the program dictates use of second as the unit for time. 

This sample run is for the all-movable control surface example of reference 1. 

The flutter equation is formulated with five vibration modes and the aerodynamic matrices 
have been calculated for a Mach number of 0, 6. A matched -point flutter solution is 
required, and the initial values of air density and the inverse of reduced frequency are 

1. 146797839 X 10 Ib-sec jj. ( se a-level density from ref. 2) and 6.5, respectively. A 

in 4 

detailed output is desired, and IPRT * 2 is input. The namelist input for a sample run 
follows. 


NAMELIST INPUT FOR SAMPLE RUN 


program no LANGLEY RESEARCH CENTER date 

CODED BY FORTRAN ■ DATA CODING FORM PAGE 0 F 

DIVISION SECTION JOB ORDER TASK NO 


■TATKM KMT 
NUMPtn 

f 

-coMTiNu.nc-N FORTRAN STATEMENT 

IDENTIFICATION 
ANO &E.QUCNCINO 

1 1 i 1 3 


F ■ A.!. 1. 



,$ N AM 

A 


CH ip 

ERF,=i 1.0,, AX MAT, = 3 , M A,C H = 0 . ,6 , . l .T .R,tf .P .0 s.0,,1 MAT.CH = 1 , ,REFS L.D = 1 . 1,4 6 7 9 7,8 3 

9 . 



1 






■ 

1 

, . ju 

jil,T.IKl,2.. 1 0,t , . . . , , 





, ! 



, 



, i 


mum 

$,NAM 

1 


S.KCl 

) = 4 ,.10 4 9, E 0,1 ., 12 *, 0 . 0 ,1,. 01 6.E,0 ,3 , 1. 2,* ,0. . 0 ,, 6, . 2.E 0,3 , 1.2 *,0 . 0 , 5,. 5 2 1 E,0 3 , 1 2, *0 

■iiipg 





J l_L ■ ■ ■ ■ 1 ■ ■ 1 . 1 1 ■ ■ ■ t . . . i 1 ■ ■ ■ E t ■ , . ■ 1 ■ ■ ■ ■ t ■ ■ , , . 1 . . .11 1 1 





1 






, 1 


mSSSSm 




C 1 

),=,7 , .;6 ,3, 1,E -r0 ,3 , , , 1 .2, *,0, . , 0, . ,4 , . ,7.4, 1, E -,0.3. . , 1,2 , * 0. .i0 , . .2. . .4.9 E. -.0 .3, . . 1. 2.*.0 , . 0 ..7-.i5.fi 2 E - fl.4 

'mJaOvjmM 


■ 


. r 5 


wmmm 





btMi 111 mu 

Bwfimi 

PHHPPP 


P 



mm 





EHjEa 


hhhhhhw 

1 , , , , 



, ! 


hhhhhhji 


■ 


N I TjM 


mmm 

' M 

■ 


, , ! 


PHiimi 


I 




■HUP 

HHI 

_ 













■ 




. . f . . . . 

mu 

1 



,, p ... 1 ........ , . . . , i 





1 

,,1,,-,,, , , , . . , . . . , 








mu 



■MS 


a — — 11 


NASA-Langley Form 67 (MAR 69) NOTE: WRITE NUMBERS 10, LETTERS I^UGZC, SYMBOLS / 
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APPENDIX A - Continued 


Note that namelist NAM1 does require an input for air density since ITROPO = -1 in 
namelist NAMATCH. 

The program output is in two parts : (a) summary on tape 7 and (b) output file. 
The output obtained from the sample run follows. 

(a) Listing of tape 7 

DENSITY « 1.1467986-07 , RFIMIN * 6.5000 


ROOT 

NUMBER 

1 

, VELOCITY * 

6316.702 # 

RFI « 

10.000 , 

NO. 

OF 

ITERATI ONS 

REQD, *> 

3 

ROOT 

NUMBER 

2 

# VELOCITY = 

29430.954 . 

RFI = 

11.350 , 

NO. 

OF 

ITERATI ONS 

REQD, « 

1 


FOUND 2 ROOTS * RFI FOR THE NEXT ROOT PR6CICTED = 22.8007 ,1$ BEYOND THE RANGE 


ITERATION NO. 1 DENSITY = I.1468E-07 S CRT ( SEA LEVEL DENS IT Y/DENSI T Y I = 1.0000E+00 

FLUTTER VEL * 6.3167E+03 AIR SPEED = 8.Q361E+03 ( VE L-AI R SPEED » *100/ VEL a -27.2204 


OENSITY * 7 . 80613 IE- 08 , RFIMIN « 13.1467 


ROOT 

NUMBER 

1 i 

VELOCITY = 7650.904 

* RFI « 

13.150 , 

NO. OF 

ITERATI ONS 

REQD, * 

1 

ROOT 

NUMBER 

2 , 

VELOCITY » 35709.560 

# RFI « 

13.800 , 

NO. OF 

ITERATI ONS 

REQO# = 

1 


FOUND 2 ROOTS# RFI FOR THE NEXT ROOT PREDICTED * 27.64C3 ,IS BEYOND THE RANGE 


ITERATION NO. 2 DENSITY * 7.8061E-08 SQRTlSEA LEVEL DENS I TY/DENS IT Y J = 1.2121E+00 

FLUTTER VEL * 7.6509E+03 AIR SPEED * 7.6504E+03 l VEL -AI RSPEED > *100/ VEL * .0065 
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APPENDIX A — Continued 


(b) Listing of output file 
SNAMATCH 

P«=RF = 0 . 15 + 01 , 

MAXMAT = 3 , 

MACH = 0 . 6 F+ 00 , 

ItroPO = - 1 , 

IMATCH - 1 , 

RPFSLO = 0 . 1146797839 F- 06 , 

UNITL = 0 . 12 c + 02 , 

t C N 0 


$ NA Ml 


0 .4049 F* 0 1 , 

0.0, 

0.0, 

0.0, 

0.0, 

O.Ot 0.0, 0.0, 0.0, 

0.0, 

0.0, 

0. 0, 

0.0, 

0.1016F+03, 

0.0, 

0.0, 0.0, 0.0, 0.0, 

0.0, 

0.0, 

0. 0, 

0.0, 

0.0 , 

0.0, 

0.0, 

0.62 c +03, 0.0, 0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0 .Of 

0.0, 

0.0, 

0.0, 0.0, 0. 0, 0.5521F+03, 

0 o 0 , 

0. Of 

0.0, 

0.0 , 

0.0 , 

0.0, 

0.0, O.Of O.Ot 0.0, 

0.0, 

0.0, 

0. 2l65 c +0<f, 

o 

. 

o 

0.0, 

0.0, 

. 

o 

• 

o 

o 

• 

o 

*• 

o 

• 

o 

o 

• 

o 

0. 0, 

0 .0, 

0. 0, 

0.0, 

0.0, 

It I, 

It It 

I, I, I, I, I, T, T, I, 

I, T, I, 

t If 

If If 

I, I, 

I, I, 

It I , 

It I, 

It I, I, I, I, I, I, I, 

I, It I, 

f I, 

It I, 

It I, 

It I , 

I, I, 

It It 

I, I, I, I, I, I, I, 1, 

If I, If 

f I, 

It I, 

I, I, 

It I f 

1 1 I , 

I, I, 

I, It It I, It I , 


0-7631^-02, 

0. 0, 

0. 1, 

0.0, 

0.0, 

0.0, 0.0, 0.0, O.Of 

0.0, 

0.0, 

0. 0, 

0 .0, 

0 .47^ 1F-02 » 

0.0, 

O.O, 0.0, O.Of 0.0, 

0. 0, 

0.0 , 

0.0, 

0.0, 

0.0, 

0.0, 

O.Ot 

0 . 249F- 02 , 0.0, 0.0, 

0.0 , 

0.0, 

0. 0, 

0.0, 

0.0, 

0.0, 

O.Of 

0.0, 0.0, 3.0, 9.7562F-03, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0 .0, 

0.0, 0.0, 0.0, O.Of 

0.0, 

0.0, 

0. 1933 c -02f 

0. 0, 

o 

o 

0.0, 

0.0, 0.0, 0.0, O.Ot 

0.0, 


0 « 0 » 0 . 0 , 0 . 0 , 0 . 0 , t , I, I, I, I, I , I, If If If If I, If If If 
1 » *» 1 » I* It It If I, I, I, I, If If I, I, If If Tf If I, If If I, 
It If If If If I, If I, If If If, If If If If If If If If I, If If If 
It If If If I, If If If If If If If If If If If If If 


LSTIFF 

= 

I , 

ss 

= 

0.16F+02 

B» 

= 

0.65 c +01 

NM 


5, 

RFIL 


0.1E+0 1, 

RFIR 

= 

0.2F+02, 
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APPENDIX A - Continued 


DEL 


0.5E-01. 

NPQOT 

= 

3 » 

NTT MAX 

= 

5, 

ND 

= 

1» 

RHQ 

= 

0.1146797839E-06 

R F l M I N 

= 

0 . 6 5 C + 01 , 

I PR T 

= 

2, 

I OPT 

- 

Ot 


$ END 



DENSITV = l. 

14679SE-0? , 

RFIMIN = 

6.5000 






E IGENVALUE S 

1.814F-03 

-1. 174E-04 

7.226 c -05 

-2.57BF-06 

4.607^-06 

-1.013R-07 

2.170E-06 

-2 .38BE-08 

9. 0R3E-07 

EIGENVECTORS 

9.941 5-Ot 
1. 75TE-02 
1. 530E-03 
-2.434E-04 
-3. 290E-04 

1 . 063E-01 
-2. 1745-03 
3. 048E-04 
-K799F-05 
-5. 137E-05 

-2. 113 F- 02 
9 . 979 F -0 J. 
-3.409F-02 
6.669E-03 
7.X55F-03 

— 3. 943F-02 
2. 9 85^-02 
-2.620E-03 
1.006E-03 
2.551 f-04 

-1. 85 ?f_02 
4.0525-03 
1.179F-02 
-4. 2^^~03 
-5.493 r -0 3 

7. 7 82 c -02 
2. 246 E- 31 
9.6575-01 
1 . 020 C -0I 
-3-1695-32 

-1 .012F-01 
-5 . 166F-02 
2. 288F-01 
-4.8 99 E-01 
1.3T4E-02 

1 .203F-0! 
B.887F-02 
-3.272F-01 
7 .513F-01 
-2.593F-02 

5.704F-03 
9.1 86E- 02 
-1.442E-02 
-3.822E-04 
-5.253F-01 

ASSOCIATED EIGENVECTORS 

1. 3595*02 -5.660EHN) 

1 . 247E-01 2.024F-01 

-5.730E-02 -1. 883E-02 

-2.049E-01 -2.8195-02 

3.544E-03 1 . 960E- 03 

6 .1 26 F*01 
1 .355F*02 
-4.6005*00 
-5.099E*00 
i .oi7 p *oo 

1.476F*00 
8. 425 c -01 
2 . 209 c -03 
4. 345E-02 
8.877F-02 

' -7.74 6 e *0Q 

4. 253 c *00 
1 . 049 c +0 1 
1.0905*01 

5. 029E-01 

4 .719 c *01 
-6.5 95 5 * 01 
-3.436 c *02 
-1.698E*02 
3.1415*00 

-1. 553F* 01 
2 .2095*01 
4.0565*01 
-4. 811F*02 
-3.277F-01 

-2.2265*01 
4.470F*01 
7.124 p *01 
- 7.4975* Q 2 
-2.077F-01 

-2.1 64F+01 
3.338 p *01 
-3.3 90E *01 
-5.459E+01 
-2.5355*02 


DERIVATIVES OF INVERSE OF FREQ. SQUARED, NUMP P R 1 
FIRST OERIVATIVF -2. 384F-05 -1.R15E-05 

SFCOND DERIVATIVE 


-3.889F-06 -2. 21 4E-07= 7.335F-06 5.585F-06* 

-1 .6475-07 

-1 .5005-07* 

-1 .106E-05 

-5.656P-06 

DERIVATIVES OF INVERSE OF FREQ. SQUARFQ, NUMBER 2 
FIRST DERIVATIVE 8.361E-06 -8.565 p -08 

SECOND DERIVATIVE 

1.503F-06 2. 10 7F-07= -2.573 c -06 2.635F-08* 

1.925E-07 

I .493F-07+ 

3 . 883F-06 

3. 5 10F-08 

DERIVATIVES OF INVERSE OF FREO. SQUARED, NUMBFR 3 
FIRST DERIVATIVE 1.682E-07 -1.909F-O9 

SECOND DERIVATIVE 

3.637E-08 6.82 8E-09= -5.175S-08 6.119'=-10 + 

9. 853 c - 09 

4.4Bl p -09* 

7. 826F-08 

I.735E-09 

DERIVATIVES OF INVERSE OF FREQ. SQUARED, NUMBER 4 
FIRST OERIVATIVF 2.029F-07 -7.093F-09 

SECOND DERIVATIVE 

-1.305E-09 -3. 5905-09= -6.243F-08 2.182F-09* 

-3.242F-08 

-3.042E-09* 

9.3555-08 

-2.731F-09 

DERIVATIVES OF INVERSE OF FR FQ . SQUARED, NUMBER 5 
FIRST DERIVATIVE -9. 752E-09 -7.242E-09 

S c CONO DERIVATIVE 

-6.749E-09 -8.536E-10= 3.001E-09 2.228F-09* 

-5.247 c -09 

-6.799F-10* 

-4.503F-09 

-2.402F-09 


SUBROUTINE LEGROOT 
DAMPING 

FIRST DERIV 

SECOND D C RI V 

RADICAL 

IN L r GUFRR c 

eigenvalue no. 

PROJECTED CROSSING 

-1. 1006E-02 

-2.2397E-03 

-1.2423E-03 


-8.6567 c -06 

4 

1.00005*03 

-2. 197BE-02 

3.7063F-04 

1 .62035-03 


3. 5923C-05 

3 

1.0167F*0l 

— 3.56 73E-02 

2.94245-03 

2. 976 9 p -Q 3 


1.1 465F-04 

2 

9.8287R+00 

-4. 3937E-02 

-8. 4452 F- 03 

-1 .44765-03 


7.7175F-06 

5 

2.2316 p *0t 

-6.4711E-02 

-I.0859E-02 

-5.4637E-04 


8. 25&4F-05 

1 

l.3622 p *Dl 

ITERATION 1 
6.5000 

9.8287 10. 

1669 13.6217 

22.3160 

1000.0003 




3.991 F-08 


6.1 08 c -03 
1.200E-01 
3.31 8F-02 
5.495F-04 
8.365=-01 


4.550FOI 
6.238F+01 
4. 1 45F*Q1 
7.655E+01 
4«463 F *02 
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EIGENVALUES 


1.7115-03 

-1.798E-04 

1 .0935-04 

-1.2455-06 

5.5245-06 

-5.6025-08 

2.6985-06 

-7.B86E-08 

8.3565-07 

-6 . 9 99 C -QR 

EIGENVECTORS 

9.758E-01 
-4. 246E-02 
3. 682E-03 
-5.B92E-04 

-7. 939E-04 

2. 142E-01 
-1. 257E-02 
1. 2805-03 
- 1.564E-04 
-2. 520E-Q4 

-6.Q39E-02 
9.930E-0I 
-5.3865-02 
1 .0505-02 
1 .145 5-02 

-6.4 09E-02 
-5.5615-02 
-2. 6495-04 
4. 330F-04 
-2.518F-04 

-1.923E-02 
1. 2845-02 
7. 4855-02 
7.344 5-03 
-7 .3005-03 

2 . 389 c -0 1 
4. 6005-01 
8. 053 c -01 
2.7375-01 
-4.0865-02 

2.736F-02 
1 .566C-02 
-9.5965-02 
1.199^-01 
5. 6365-03 

1.8735-01 
3.245^-02 
-7.5555-01 
6 .070 c - 01 
-6. 979*=-03 

-4.1705-03 
1.1735-01 
5 . 776F-02 
-3.4365-02 
-3.6875-01 

-2. 7165-03 
-2.201 c -0! 
-6.003F-02 
5.2705-02 
8.899R-0I 

ASSOCIATED EIGENVECTORS 

1 . A3 25*02 -1. 554E+01 

3.467E-01 3. 427E-01 

-1.4B4E-01 -2.961E-02 

-5. 2B7E-01 -3.823F-02 

9.317E-03 3. 298E-03 

1.012F+02 
8. 8905+01 
-5 .5045*00 
-6.799F+00 
1.0795+00 

1.4645+01 
6, 152E+00 
-6. 233E-0I 
-1.2075+00 
1.430 c -01 

- 1. 070 c +01 
5.314 c +00 
1.5775+01 
4. 578 c +01 
9.543 c -01 

4.980^+01 
-6.276F+01 
-2.5105+02 
-3. 4405+02 
-9.0245-01 

1.3605+01 
-3.2735+01 
-5.560 c +01 
1 .4275+02 
-8.4355+00 

-6.253 c +01 
1 .0965+02 
2.195F+02 
-7.609^+02 
1.7865+01 

-1.807F+01 
2. 8695 + 01 
-3. 578F+01 
-2. 6195+01 
-1. 7405+02 

-8.2095+01 
1. 1885+07 
-6.7385+01 
-1.2475+02 
-5.431C+02 


DERIVATIVES OF INVERSE OF FREQ. SQUARED, NUMBER 1 
FIRST DERIVATIVE -3.745E-05 -1.930E-05 

S ECONO DERIVATIVE 


-4.298E-06 -5.0975-07= 7. 6055-06 3. 919^-06 + 

-4. 5605-07 

-3. 804F-07+ 

-1.145F-05 

-4.0495-06 

DERIVATIVES OF INVERSE OF FREQ. SQUARED, NUMBER 2 
FIRST DERIVATIVE 1.3995-05 1.0425-06 

SECOND DERIVATIVE 

1.907E-06 5 . 044E-0?= -2.B41E-06 -2.116E-07+ 

4.67B c -07 

3.775F-07+ 

4.280F-06 

3.3855-07 

DERIVATIVES OF INVERSE OF FREQ. SQUARED, NUMBER 3 
FIRST DEft IVATIVE 4»22B c -07 3.084F-08 

SECOND DERIVATIVE 

l .0905-07 9.4985-09= -8.585E-08 -6.2615-09+ 

6.5705-08 

5.722F-09+ 

1.291 c -07 

1.0045-0B 

DERIVATIVES OF INVFRSE OF FREQ. SQUARFD, NUMBER 4 
FIRST DERIVATIVE 7.3245-08 -2. 6035-08 

SECOND DERIVATIVE 

-6 . 581E-OB -3. 57 6E-09= -1.487F-08 5.286'=-09 + 

-7.3165-08 

-1.743F-09+ 

2.222F-06 

-7.1 19F-09 

DERIVATIVES OF INVERSE. DF FREQ. SQUARED, NUMBER 5 
FIRST 05RIVATIV P -3.4195-08 -1.1035-08 

SECOND DERIVATIVE 

-7.834E-09 -1.4565-09= 6.943E-09 2. 2405-09+ 

-4. 32 85-09 

-l .049E-09+ 

-1 .0455-08 

-2. 647^-09 


SUBROUTINE LEGROOT 
OAHPING 

FIRST DERIV 

S c COND DERIV 

RADICAL 

IN IFGUERRF 

EIGENVALJF NO. 

PROJECTED CROSSING 

- 1. 01 41E-02 

6. 3583 5-03 

9 .46215—04 


5.0024 e -05 

3 

1 . 1 284 c +01 

-1 . 13 90E-0Z 

1.09895-02 

1.9996F-03 


1.4354F-04 

2 

1.08015+01 

-2. 92335-02 

-8. 85695-03 

-1.5577E-03 


3.2907 c -05 

4 

1.49465+01 

-8. 376 OE -02 

-1. 6632 c -02 

-3. 08885-03 

- 

■4. 9108F-0 5 

5 

1.0000^+03 

-1.05095-01 

-1.3580 E -02 

-1. 1563E-03 


6.2920E-05 

1 

2.3098P+01 

ITERATION 2 
9. 8500 

10.8007 11. 

2838 14.9460 

23 .0979 

1000.0000 




to 

w 
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to 


EIGENVALUES 


I. 6746-03 

-1. 9B4E- OA 

1.235E-0A 

-6. 779F-09 

5.976E-06 

-2.273F-38 

2.737f-06 

-1 .0476- 07 

7.9956-07 

-8.1176-08 

E IGENVECTDRS 

9.7186-01 
-5.209E-02 
A. 516F-03 
—7 . 23 A F-04 
-9. 7396—04 

2. 2916-01 
-1.729E-02 
1.71AE-03 
-2. 206F-0A 
-3. A22E-0A 

-6.531E-02 
9.9 08 E-01 
-S.76A=-02 
1.117E-02 
1.233E-02 

-8.207E-02 
6.105E-02 
-7.1586-03 
1. 710E-03 
1 .2256-03 

2.867E-31 
A.926F-0I 
6.596F-01 
2 . 909E-01 
-3.507F-02 

-8.23ZF-32 
-1.877F-31 
-3.135E-01 
-l . 1 61 E-Ol 
1.986E-02 

1 .605 e-01 
-5.950E-03 
—7 .41 76—01 
A. 7986-01 
-5.3226-03 

-8.086E-02 
-9.165E-03 
3 .4436-01 
-2. 6206-01 
-7.259F-03 

-7.9116-03 
2. 4116-01 
1 . A 01 F-Q 1 
-8 • 150E-02 
-8.000E-01 

-3. 6736-03 
-1.266F-01 
-3.8836-02 
3.31BF-02 
5.069F-01 

ASSOCIATED EIGENVECTORS 

1.4646*02 - 1. 59BE* 01 

4. 3496-01 A. 015E-01 

-1.860E-01 -3.6A9F-02 

-6.634 c -01 -5. 222E-02 

1.171E-02 3.9996-03 

1 .125 c *02 
7. 873E»0! 
-5.6986*00 
-9.779F*00 
1 .0506*00 

3.79ae*00 
-A. 6226*00 
A. 8D3E-02 
-2.072F-01 
-2. 3666-03 

-4.074 c *0l 
5.1 09R+01 
2 . 0 A 6 F*O 2 
3.3366*02 
1.098 c *00 

-2.557F*01 
2.2386*01 
8.561 F*0! 
1.663E+02 
1.2 91 c *00 

6.6626*01 
-1. 1346*02 
-2.3556*02 
6.6076*02 
-2.3456*01 

3.8226*01 
- 7 • 557 e * 01 
-1 .4376*02 
3.5096*02 
-2.1156*01 

- 6 . 9?3F*01 
1 . D3 46 *02 
-7.2846*01 
-9.8506*01 
-4.7676*02 

-6.6256*01 
9.5906*0! 
-3.6116*01 
-1.016F*02 
-3 . 73 2F *02 


0ER1VATIVFS 06 INVERSE DF FR^q. SOUARED, NUMBER 1 
FIRST DERIVATIVE -4.162 c -05 -1.986 c -05 

SECOND DERIVATIVE 


-4.4706-06 -6.6316-07= 7.7076-06 3.677F-06* 

-5.8666-07 

-5 .033 F-07 * 

-l .1 596-05 

-3.837F-06 

DERIVATIVES DF INVFRS 6 OF FRFQ. SQUARE 0 . NUMB 6 R 2 
FIRST DERIVATIVE 1.5886-05 1.589 c -06 

SFCDND DERIVATIVE 

2.0826-06 6. 554E-07* -2.9416-06 -2.943F-07* 

5.92 76-07 

4.993F-07* 

4.4316-06 

4.504F-07 

DERIVATIVES OF INVERSE OF FReq. SQUARED. NUMBER 3 
FIRST DERIVATIVE 5.287F-07 3.896 C -Q8 

S ECONO derivative 

1.1256-07 7. 736 c - 09 = -9.790E-08 -7.214E-09* 

6 .3206-08 

3.8786-09* 

1 .4726-07 

1.1O7E-08 

0FR1VATIVFS OF INVERSE DF FREQ. SQUAREO, NUMBER 4 
6 IRST DERIVATIVE l.D90 c -08 -2.7816-08 

SECOND DERIVATIVE 

-6.387E-08 -2.579E-I3- -2.0136-09 5.150E-09* 

-6.4826-08 

1 .420 6-09* 

2.969F-09 

-6. 8286-09 

DERIVATIVES (IF INVERSE DF FREQ. SQUARED, NUMBER 5 
ETRST DERIVATIVE -4. 1H8E-0B -1.256E-08 

SFCONO DERIVATIVE 

-8.3796-09 — 1 . 7726-09= 7.755E-09 2.325 e -09* 

-4. 463E-09 

-1.2986-09* 

-l .16TF-08 

-2.8006-09 


SUBROUTINE LFGRCOT 
OAMP ing 

FIRST 06 R IV 

SECOND deriv 

RADICAL 

IN LEGU 6 RRC 

C TGFNVALUF NO* 

PPDJFCTED CROSSING 

-5.4B84E-05 

1. 28 74 c - 02 

1.99626-03 


1 .65866-34 

2 

1.08046*01 

-3. 80436-03 

6. 8553 c -0? 

1.53216-04 


4. 7579 c -05 

3 

1.13526*0! 

-3.32426-02 

-1.00066-02 

-9.06776-04 


6 . 5447 c -05 

4 

1.55276*01 

-1.01536-01 

-2. !024=-02 

-5.48346-03 


-1 . 14 72 6 — D4 

5 

1.00006*03 

-1.18556-01 

- 1. 4812E-02 

-1.44956-03 


4.7553 E-0 5 

1 

2.79916*01 

ITERATION 3 
10.8000 

10.8043 1 

1.3515 15.5271 

27.9913 

1000.0000 




FLUTTER EIGENVALUE NO. » 2 t eigenvalues 

1.6737F-03 -1.9841E-04 1.23516-04 -6.77866-09 5.97596-06 -2.2734 c -08 ?.7375=-06 -I.0469E-0? T. 99496-07 -8.T170E-08 


RF1 FOR PREOTCTFO CROSSINGS COPP C SPONO TO EIGENVALUES NUMBERS 
2 3 A 1 5 


, RFI 


ROOT NUMOFR 


1 , VELOCITY 


6316.702 


10. 800 


NO. OF ITERATIONS R c q 0 
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E IGFNVALUF S 


1. 6506-03 

-2.0946-04 

1.326F-04 

9.7186-07 

6. 284=-06 

-1.6796-1 0 

2.734F-06 

-1 .1996-07 

7.7526-07 

-8.836F-08 

6 IGFNVFCT ORS 

9.657F-01 
-5. 795E-02 
5. 0196-03 
-8.0546-04 
-1.083=-03 

2. 5216-01 
-2. 1486-02 
2.0926-03 
-2.7816-04 
-4. 216F-04 

-8. 3355-02 
9.899 c -01 
-6.0265-02 
1.1706-02 
1.2806-02 

-8.0375-02 
-5.2 B2 6-02 
-5.4556-04 
3. R69F-04 
-1. 8506-04 

3.0686-01 

5.3496-01 

6.6406-0! 

3.1536-0! 

-3.8356-02 

1 .2496-Jl 

1.736F-01 
1. 5316-01 
9.3806-02 
-6. 2096-03 

1.6796-01 

-1.7146-02 

-8.0606-01 

4.9836-01 

2.6486-03 

4.3386-02 
-1 . 5756- 02 
-2.4076-01 
1.1576-01 
-9.8176-03 

-l. 1Z5F-02 

2.7616-01 
1.8726-01 
-l .0256-01 
-8.7996-01 

-5.3606-03 
-T. 3656-02 
-1.3 56 F-02 
1.7265-02 
3.1306-01 

AS50CIATFD 6 IGENVEC TOP. S 

1.4836*02 - 1. 835F* 01 

5.0006-01 4. 3386-01 

-2. 1166-01 -3. 844E-02 

-7.5396-01 -5.2326-02 

1.3366-02 4.3035-03 

1 .1 706*02 
7.3406*01 
-5 .6755*00 
-1.0175*01 
1 .0166*00 

1.8036*01 
3.6526*00 
-5. 9436-01 
-1.3896*00 
1.0446-01 

-4.6326*01 
4.9866*01 
1.9856*02 
3. 7295*02 
1.7966*00 

5.0426*00 

-1.3716*0! 

-5.7656*01 

-8.3276*01 

1.9036-01 

8.2156*01 
-1.4706*02 
— 3.989F*02 
7.1 346*02 
-4.0356*01 

— 1 . 2255*01 
2.1245*01 
5.4846*01 
-1.7396*02 
-6.9756- 02 

-8. 8436*01 
1.3196*02 
-8.1966*01 
-1 .249F *02 
-5.68BF+02 

-5.4736*01 

7.9516*01 

-1.8996*01 

-8.6676*0! 

-2.767F*02 


D6RIVATIVFS OF INVFRS r OF FREQ. SQUARED, NUMBER 1 
FIRST OFPJVATIVF -4.4115-05 -2.025F-Q5 

S=CDND OF R f VAT I VF 


-4. 586 c -0fc -7.7676-07= T.77zc_ 06 3.5696-06* 

-6. 755=-07 

-5.9536-07* 

-1 .168F-05 

-3.750F-06 

DERIVATIVES OF INVERSE OF FRFQ. SQUARED, NUMB 6 R 2 
FtRST DERIVATIVE 1. 706 c -05 3.979F-06 

S c COND DERIVATIVE 

2.2016-06 7.6715-07= -3. 0066-06 -3.488 c -07* 

6. 7826-07 

5.9066-07* 

4.5296-06 

5.2 53 6-07 

DERIVATIVES OF INVERSE DF FREQ. SOU AR p O* NUMBER 3 
FIRST DERIVATIVF 5.904 c -07 4.3066-08 

SFCOND DF R IVATI V5 

1.117F-07 7.2746-09= -1.040E-07 -7.5886-09* 

5.9326-08 

3.477F-09* 

1.5646-07 

1.1436-08 

DERIVATIVES DF INV C RSF DF FR6Q. SQ1J6R50, NUM0FR 4 
FIRST DERIVATIVE -2.3095-08 -2.7525-08 

SECOND DER ! VAT 1V5 

-5.9476-08 1.2226-09= 4.0686-09 4.8506-09* 

-5.7396-08 

2.7176-09* 

-6.1446-09 

-6.3456-09 

DFRIVATIvFS OF INVERSE 06 FR 6 0 . SOUARFD, NUMBER 5 
FIRST DERIVATIVE -4.660F-08 -I .3606-08 

SFCOND DERIVATIVF 

-8.8276-09 -2.0206-09= 8.2126-09 2.39fe6.(59» 

-4.6816-09 ' 

-1.5006-09* 

-1 .236E-08 

-2. 9166-09 


SUBROUTINE LFGROOT 
DAMPING 
7.33086-03 
-2.67275-05 
-4.38706-02 
-1 . 1 ? 98E-0 1 
- 1 . 26935-01 


FIRST OF R IV 
1, 39695-02 
6. S553=-03 
—1.04386-02 
-2.4393 c -02 
-1.56656-02 


SECOND DFRIV 
2. 06456-03 
-1 . 30136-04 
—6. 83526—05 
- 6 . 83 7 OF— 03 
-1.66096-03 


D A 0 1 C A L IN L6GU C #R6 
1.80556-09 
9.69926-05 
7.89696-05 
' -1.84266-39 
3.45926-05 


EIGENVALUE NO. 

2 


4 

5 
1 


PRDJ6CT 60 CROSSING 
1. 08046*01 
1.13546*01 
1.62876*01 
1.00006*03 
3. 29315*01 


ITERATION 1 
11.3500 


10.8044 11.3539 16.2369 32.9306 1000. OOQO 


FLUTTFR EIGENVALUE NO. = 3, 6IG=NVALU6S 

1.65016-03 -2.09446-04 t. 32566-04 9.71805-07 6. 28-366-06 -1.6794E-10 2.73405-06 -l. 19946-07 7.751 7F-07 -8.83566-08 


RF I FOR PRF01CT50 CROSSINGS COR RFS POND TO EIGENVALUES NUMB5 B S 
2 3 4 1 5 


ROOT NUMBPR 2 , VELOCITY, = 29430.954 . R F I - 11.350 . NO. OF ITERATIONS R6QD , = 1 

M 

CJl 
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P IGFNVAU/f S 


1. 370 c -02 

-3. 2698-04 

2 .49R c -04 

2.787F-05 

1.35 4 c -05 

3.19 ,c -07 

2 ♦ 090 p - 06 

-2.109F-07 

4.017P-07 

-! .949 c -37 

c 5GFNV c CTORs 
-7.385F-01 
9. 177 c -02 
-7. 764F-03 
1.291F-03 
1 • 7 0? c -03 

-6. 5448-01 
1. 3348-01 
187E-02 
1 . 8 49 F- 03 
2. 515F- 03 

1 .431F-01 
-9.352 C -Q1 
6 .927 c -0? 
-1.318=-02 
-1.801F-02 

2. 277F-01 
- 2 . 1 80 c -01 
2 .0 70 r -0? 
-3. 750 c ’-')3 
-4. 202F-03 

5. ?62 r -01 
7.777 r — 03 

3.098<-_ 0 i 

2. 977 c -01 
-1 . 7?8 c -02 

R.4 e 9 c -02 
1 .06RF-01 
I .4 56 c -03 
2.?96 r -02 

1 .44TF-03 

-7 . 224 p -02 
6 . T 99 C — 02 
5. 167F-01 
-1 .S84F-01 
3.491 F-02 

1 .166F-01 
—9 . 89R C - 02 
-7.577F-01 
2 *982 r — 0! 
-3 .929F-03 

4.420F— 02 
- 1 . 688 F -01 
-2.R4IF-01 
1 . 3 74 e- 01 
?.591 c -0! 

- 5. 06 IF- 02 
2 . 886 F- 0 ! 
5-. 381 e -01 
-2.085F-0' 
-5. 632 C -0I 

ASSOC [ATFp e IGENVFCrOBS 

- 1 • 64 0 p * 02 7. 92! c * 01 

-1.7B6 P *00 -6. 441F-01 

6.4I6F-01 -5. 12 RF -02 

2 . 223 p *00 -2. 724 c - 01 

-4.304 p -02 -1. 884F-03 

-1.736P+02 
-3 . 8 32 c * 0 ! 
5 ,363 p *00 
1.36 OF *01 
-7.013F-01 

9.738F*30 
1 . 201 F *0 1 
-1. 144F + 00 
-2. 2 ? 5 r * 00 
2.356 c — 01 

-4.635^*0* 
3.394 p *01 
l . 25 3 C *02 
4.01 TF*02 

1 .?0R P *D9 

-6.132 p *30 
-2.9605*00 
-1 . 566F*ol 
-3.146 p *0 1 
9.164F-02 

-7.11fcF*01 
1.14JFFQ2 
2 . 76 7F » 02 
-3 . R 02 p *02 
8 .601 r *0l 

-l .0524*02 
2 ♦ 003 c * 02 
5.006F+02 
-6 .008 C *Q2 
2 . 3 92.F* 02 

— 4. 6 1 8 P ♦ 01 
8.168 p *01 
1 . 210 F +02 
-1.364? *02 
l.787F*0! 

3. F68 c *02 
-5.696F*02 
6 . 060F*0 1 
5. 5 47F* 07 
1. 622 c *0? 


OFRrVftTIVFS OP 1 NVFRS C OP F«FQ. SQUARED, MUMBFP 1 
FIRST n c R!V 4 TrV r - 7 . 06 Z c -0 5 -2 . 984 p -05 

S F CONn Dfr IVflTIVE 


-6.194F-06 -4. 19 5F-06= 8.66 cc -06 3.661 r -0fc* 

-1.764F-06 

-3 . 568 P-06 * 

-1 . 31 O c -05 

— 4.2 8 9 C — 06 

DFBIVATTVFS OF !NVFPS C OF FRFQ. SOUAP r D, NUMBER 2 
FIRST DFR 1 VAT [ v r 3.169F-35 i. 150 c -05 

S c CflND D e P IVATC VF 

3. 75 8 p -06 4.18 ? p -06 = -3 .809^-06 -1.4UF-06* 

1 .736 c -06 

3.557F-06* 

5.911 c -0fc 

2.078F-06 

DERIVATIVES OF INVFRSF OF FREq. SQUARED, NUMR c R 3 
FIRST D C R!VATIV” 1.132F-06 9.20 ie-oh 

s c rnMo derivative 

1.147F-07 1 . 41 OF- 08 = -1.389F-07 -!.129=-08* 

4 . 494 n - 03 

9.646 c -09* 

2.3R6F-07 

1.575F-08 

Dfr IVAT IVFS df INVFRSE 9F FR c Q. S3IIAPFD, NUMR c p 4 
F rRST D?P I VAT1 VF -1. 99‘3 c -07 -3.R85F-09 

S c CONO D c R I VAT I v c 

- 1 . 28?F -08 R . 12 3F-37= 2.442 C -0P 4.746 r -!0* 

-6. 379F- 10 

7.56Z c -09* 

-3 . 660F-0B 

8.480F-11 


derivatives pf inverse of f»fq. soiiarfd, ku^rpr 5 
FIRST DEP I V AT ! V c - 1 . 160^-07 -?. 50 ?f -08 

S r rOND OCR I VAT! VF 

- 2 . 378 F -08 - 7 , 645 E -09 = 1 . 423 F-OR 4 . 297 c - 09 * - 1 . 655 c -C 8 - 6 . 230 F- 09 * - 2 . 146 F -08 - 5 . 712 p -09 


SUBRDUTINF LPGROOT 
DAMPING 
*.U 56 f- 0 l 
3 . 0272 R -02 
- 1 . 0090 F -01 
- 2.78 7 CF -01 
- 4 . 90 ! 9 P- 0 ! 


FIRST OF R ! V 
3 . 1984 r -02 
5.4797F-02 
l467 c -02 
. 4091 e-O 1 
- 2 . 28 74 p- 0 ! 


5 FC PNO D c P l V RADICAL 

6 . 9776 F -03 
- 1 . 6 R 06 R -04 
I. 0844 F -03 
- 7 . 6580 F- 0 ’ 

-1 ■ 80 l 7 C — 0 1 


IN L C GU C RR C c ir.PNVALU c 
2 . 3 R 16 F -04 
3 . C 1 14 c -05 
3 .4090 p -04 

6.657te-9^ 

3 . 5996 F -02 


NO. PRHJFCTFT CROSSING 

2 9 . 07 ! O c *00 

3 1 . 1 ! 9 X c *01 

4 ?. ? 801 c * 0 ! 

1 l . 0 300 c *03 

5 1 .0000 r *03 


iteration 1 

16.3000 


9 . 0 TJ 0 U .1915 22.8007 1000.0000 1000. 0000 


FOUND 2 POOTS, P F I FOR TH c NEXT ROOT ppcptCTCD 


22.800 


IS RFyfWD THF R 4 NGF 
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ITERATION NO. 
FLUTTER V5L = 


1 DENSITY = 1.14685-07 SQRT(S'=A LEVEL 05NS ! T y /O c NST TV ) = 1.0000^*00 

6- 3167E *-03 AIR SPF*D » 8.03615+03 I VFL-A! R SPEED) *100/VFI » -27.2204 


R50UC6D FREQ. = 10.803 

••The DERIVATIVES are H.R.T. SQRTtSFA LFVFl DENS! TY/OENSITY) •• 

OFRIV. OF REDUCED FREQ, = -9.48D0F-02 

OERIV. OF FREQ**2 = -4.T296E+02 

DERTV. OF VELOCITY - 6.2823F+03 

SECOND OERIV. OF RF - 1.9344E-01 , SECONO OERIV. OF FREQ»*2 = 


1 « 01 75E * 03 , SECONP OERIV. OF VEL= 6.0119^*01 


DENSITY = 7.8061 31E-0 8 , RFIMIN = 13.1467 


1.669E— 03 

-1.643E-04 

1 .248E-04 

-9. 9955-10 

6.0075-06 

-1.774F-08 

2.7355—06 

-8.727E-0B 

7.9596-07 

-6. 7S2E-08 

EIGENVECTORS 










9.757E-01 
-5. 355E-02 
4.662F-03 
—7 .412E- 04 
-1.9035-03 

2. 120E-01 
- 1. 5 82 E-02 
1.555E-03 
-2. 032E- 04 
-3.1195-04 

-T. 3325-02 
9 .930E-01 
-5 .8335-02 
1.I27E-02 
1.246E-02 

-6. 663F-02 
2.038F-02 
-4.163E-03 
1 . 0705-03 
6.3 82 E-04 

2.9505-01 
S. 126F-01 
6.8555-01 
3.0215-01 
-3. 7165-02 

-5.938F-32 
-1 .4075-01 
-2.377F-01 
-8. 7125-02 
1.539E-02 

1.687F-01 
-5.7445-03 
-7. 8055-01 
5.0695-01 
-4.0545-03 

-5.9655-02 
-7.547F-03 
2 .5 27E-01 
-1 .9465-01 
-6.4195-03 

-4.900E-03 
2.692E-01 
1.39BE-01 
-B.650E-02 
—9. 4 02F-01 

-6.7675-03 
4.B42E-02 
5. 1185-02 
-2.275F-02 
-1.046 E-Ot 

ASSOCIATED EIGENVECTORS 









1. 472E+02 
4.54BF-01 
-1. 902E-01 
-6. 7715-01 
1. 203E-02 

-1. 651E+01 
3.260E-01 
-2.646E-02 
-2.924E-02 
3. 092E-03 

1.131E+02 
7.799F+0I 
-5 .685E+00 
-9.808=+00 
1.041 c +00 

6.579E+00 
-1.407E+00 
-l .3285-01 
-4. 6 7 2F-01 
2.941 E-02 

-4. 357E+01 
5. 2935+01 
2.1075+02 
3.513 c +02 
1.311E+00 

-I ■ 9 99E+01 
1.679F+01 
6.3 80F + 0 1 
1 .275F+02 
1.052E+00 

7. 1925+01 
— 1. 245E+02 
-2.558E+02 
7 .0065*02 
-2. 7455 + 01 

2.9035+01 
-5 . 839E+01 
-1 . 1 00 e ♦ 0 2 
2 • 611E+02 
-1.7165+01 

-9. 691 E+01 
1.429E+02 
-7. 9465+01 
-1.43ZE+02 
-6. 094E+02 

-5. 1795+00 
5 . 932F+00 
1.6036*01 
-1.216F+01 
1.6B46+01 


DERIVATIVES OF INVERSE OF FRRQ. SQUARED, NUMBER 1 
FIRST DERIVATIVE -3.463E-05 -1.356E-05 

S c COND DERIVATIVE 


-3.0735-06 -3. 8955-07+ 5.268 c -06 2.0635-06+ 

-4.3755-07 

-2.941F-07+ 

-7.9045-06 

-2.1 586-06 

DERIVATIVES OF INVERSE OF FREQ. SQUAR50, NUMBER 2 
FIRST DERIVATIVE 1.331F-05 1.121E-06 

SFCONO DERIVATIVE 

1. 4605-06 3.8265-07= -2. 0245-06 -1.7055-D7+ 

4.413F-07 

2. 9195-07+ 

3. 0436-06 

2.613E-07 

DERIVATIVES OF INVFRS5 DF FREQ. SQIJARCO, NUM8FR » 
FIRST DERIVATIVE 4.4145-0? 2.692F-08 

SfCOND DERIVATIVE 

7.6655-08 4. 30 55-09= -6.713F-08 -4.094 e -09 + 

4.2915-08 

2.120E-09+ 

! .0096-07 

6.279R-09 

06R1VATIV.ES 06 INV6RSF DF FR C Q. SQUA05D, NUMR C R 4 
FIRST DERIVATIVF 6.3806-09 -1.895E-0B 

SECONO DERIVATIVE 

-4.3105-08 -5.5735-11= -9.704F-10 2.862 c -09+ 

-4. 35 75-08 

e. 841F-I0 + 

1 .44QF-D9 

-3.8226-09 


DERIVATIVES OF INVERSE OF FREQ. SQUARED, NUMBER 5 
FIRST DERIVATIVE -3. 516E-08 -8.626F-09 

SECOND DERIVATIVE 

-5. 7765-09 -1.0135-09 = 5.347F-09 1.212E-09+ -3.089 c -09 -7.4405-10+ -8. 0345-09 ’-1.581E-09 


SUBROUTINE IFGROOT 
DAMPING 
-8.0085E-06 
-2.9531 E-03 
-3. I908E-02 
-8. 4835E-02 
-9. 8396E-02 


FIRST DER.1V 
8 .9845 F-03 
4.698l=-Q3 
-6.85 34 F-03 
-1,4585 5-02 
-l.0165 c -02 


SECONO E>FRIV 
1 . 1500 c -0 3 
6. 39735-05 
1 -4.9126E-04 

-3.1 776E-03 
-B.3629E-04 


RADICAL IN LEGUERR5 
B. 0731F-05 
2.2261 5-35 
3. 1 293F-05 
- 5. 68 47E-05 

2.1048 f-05 


EIGENVALUE NO. 

2 

3 

4 

5 
1 


PRDJECTFD crossing 
1 .315! E+01 
1 . 37769*01 
l . 8854E*01 
t. 00005+03 
3.4597E+01 


APPENDIX A — Continued 



ITERATION 1 

13. 1500 13.1 50*3 13.7759 18.H540 34.5970 1000.0000 


FLUTTER EIGENVALUE NO. = 7, C 1G=NV41UFS 

1.66945-03 -1.64266-04 1.24815-04 -9.9954=-10 6.90676-06 -1.773B6-08 2.7351F-06 -8.72725-08 7.9590F-07 -6.7520 C -0B 

RFJ FOR PRFOtCTTO CROSSINGS CORRESPOND TO FlGFNVfiLtJ c S NUMB r PS 
2 3 4 1 5 


RtUlT NUMBER 1 , VELOCITY = 7650.904 . RF! = 13.150 , NO. OF ITERATIONS l c 00» = 1 


£ IG C NV4LUFS 


1.6465-03 

-1. 73ZF-04 

1 .338E-04 

8.1305-07 

6. 310E-06 

6.459=-10 

2.7305-06 

-9.954F-08 

7.7185-07 

-7.3355-08 

c IG5NV5CTnRS 










9. 772 e-01 
-6.0125-02 
5. Z36E-03 
-8. 3245-04 
-1.1265-03 

2.0305-01 
- 1.774E-02 
1. 7335-03 
-2.2935-04 
-3.487F- 04 

-8 .0015-02 
9.906E-01 
—6 • 028 E— Q2 
1 .1625-02 
1.2915-0? 

-7.497«-02 
5. 135E-02 
-6 .2225-03 
1 .433E-03 
1. 0825-03 

3. 153E-01 
5.4565-01 
6.65QE-01 
3.194E-01 
-3. 93 8 c -02 

1.067E-01 
1.4895-01 
1.3095-01 
fi. 0305-02 

-5 •? 99 e-0 3 

1.674E-01 
-1.8735-02 
— S.094F— 01 
4.9575-01 
1 .8805-03 

4.2 96F— 02 
-1 .42QE-OZ 
-2.3465-01 
1 .1675-01 
-8.2435-03 

-7.202F-03 
2.6875-01 
1 . 665E-01 
-9.538F-02 
-8.9455-01 

-9.8335-02 

9.9675-02 

9.2855-02 

-4.3975-02 

-2.6515-01 

ASSOCIATED EIGENVECTORS 









1 .496E*02 
5. 0975-01 
-2. 1495-01 
-7.662=-0l 
1.3595-02 

- 1. 424 c * 01 
3. 6605-01 
-3. 3675-02 
-4.927E-02 
3. 6985-03 

1.189F+02 

7 . 269E * 01 

-5 .6 90F *00 
-1 . 02 B c *ai 
1 .013E*00 

3. 622E*00 
-3.9766*00 
5.434E-02 
“I .6 525-01 
— 1.152 E— 0? 

-4.626E+01 
5.0235*01 
1.9965*02 
3. T57e* 02 
1.791 F*oo 

4.5015*00 
-1.1745*01 
—4. 9085*01 
-7. 205 c *o 1 
1 .3585-01 

8.2545*01 
-1.4756*02 
-3.0015*02 
7.1305*02 
-4 • 1 01 E *01 

-1 .7945*01 
2 .4065*01 
5.8505*01 
-1 .7475*02 
1 ■ 7 29 c *00 

-1.0455*02 
1.5495*02 
— 7 ♦ 756 F *01 
-1.5375*02 
-6. 2496*02 

1. 1 75 c * 01 
-1.8905*01 
3.1395*0' 
1.0675*01 
1.2375*02 


derivatives of inverse of frf«. sduared, numbfr i 
FIRST DERIVATIVE -3. 6665-05 -1.3945-05 

SECOND DERIVATIVE 


-3. 1685-36 -4.5335-07= 5.313=-06 2.005 c -06* 

-5 .0575-07 

-3.479F-OT* 

-7 .9765-06 

-2.1115-06 

0FRIVAT1VFS OF INV5RSF OF FR60. SOUARFD. NUMBF R 2 
FIRST DERIVATIVE 1.428E-Q5 1.390 c -06 

SECOND 05 R I VATI VF 

- 1.5S15-06 4.4685-07= -2.0705-06 -2.015E-07* 

5.072F-07 

3.453F-07* 

3.113E-06 

3 .0306-07 

DROIvATIVES OF INVER5F 96 FREQ, SQUARFD. NUMBER 3 
FIRST D5RIVATIV5 4.9115-07 2.9625-08 

SFCOND DERIVATIVE 

7.6066-08 4.06 SF-07 = -7. 1176-08 -4.293E-C9* 

4.032E-C8 

’ .8975-09* 

1.069F-07 

6.4655-09 

DERIVATIVES OF [NVFRSF OF FRFQ. S0UAP5D, NUMBER 4 
FIRST DERIVATIVE -2.0T2 c -os -1. 8715-08 

SECOND 05 R IVATtvE 

-4.012E-0B 7.4006-13= 3.0025-09 2.7125-09* 

-3.860E-08 

1.5755-09* 

-4.517E-09 

-3.547F-09 


r»FR IVAT IVES OF INVERSE OF FREQ. SQUARED, NUMBER 5 
FIRST OEBIVATIVE -3 . 901 e-OB -9.328F-09 

SECOND 0FR1VATIVE 

-6.093F.Q9 -1.15 15-03= 5. 6545-09 1.3526-09* -3.249E-09 -8.5695-10* -R.497E-09 -1.6465-09 


SUBROUTINE LFGROOT 
DAMPING 
6.07 7 3 C — 03 
1. 0237E-04 
— 3 . 64 5 7 = -0Z 
-9.50375-02 
-1.05195-01 


FIRST DFRIV 
9. 7437E-03 
4.6870=-03 
-7, 1 307 c - 03 
-1.6BB9E-02 
-1. 0747=-02 


S c CONO 05RTV 
1 .1 883E-03 
-0.60305-05 
-3. 7288E-04 
-3.94945-03 
-9.5653F-04 


RADICAL IN LFGUFRRE 
8.77 14F-05 
2.19775-05 
3.7253F-05 
-9.0099F-05 
1 . 48855-05 


EIGENVALUE NO. 

2 

3 

4 

5 
1 


PROJECTED CROSSING 
1.31516*01 
1 . 37785 *01 
1 .9773 c *01 
I. 00005*03 
4.10645*01 
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ITERATION l 
13.8000 


11.1511 


13.7782 


19.7731 


41.0639 1000.0000 


FLUTTER EIGFNVALUF NO. = 3, *1 GE NV ALUES 

1 , 64b25 -03 —1.731b E-04 1.33775-04 8.1299^-07 6.30985-06 6.4591^-10 2.73035-06 -9.9540E-08 7.7181F-07 -7.3350=-08 


RFI F 08 PPE01CTE.0 CROSSINGS CORRESPOND TO EIGENVALUES NUMBERS 
2 3 4 1 5 


POCT NUMBER 2 , VELOCITY = 35709.560 , RF! = 

13.800 , NO. 

OF ITERATIONS 

REQO, - 1 





6 I GFNVALU* S 

1.3656-03 -Z. 7026-04 2.533F-04 

2.375*-05 

l.055*-05 

2,6?0*-07 

2.0916-06 

-1 .736F-07 

3. 96 8 F -07 

-1.626F-07 

* IG6NV6CT0RS 

-8.583*-01 -4.6586-01 1.715F-01 

l . 2396-01 1.0826-01 -9.418*-01 

-1.0666-02 -9. 6666-03 7.019F-02 

1.7346-03 1.4996-03 -1.327*-0Z 

2. 310*-03 2.0436-03 -1.5196-02 

1 .998 *-01 
-1.951*- 01 
1.829F-02 
-3.2996-03 
-3. 726*-03 

5.171*-01 
7. 295 r -0l 
3.094 *—01 
2.977*-Q! 
- 1. 72 4F-02 

7.622*-02 
9.7146-02 
4. 9796-03 
2.265*-02 
9. 902 E-04 

— 1 . 339F- 01 
1.156F-01 
8 .873F-01 
-3.438F-01 
1.7176-02 

-2.992F-02 

3.068F-02 

2.323E-01 

-7.952*-02 

2.7186-02 

&.131E-02 
-3.2356-01 
-6. 2436-01 
2.3S2F-01 
6.1086-01 

2. 691 *-02 
-8.5396-02 
-2. 1576-01 
7.416E-02 
1.0946-01 

ASSOCIATED EIGENVECTORS 

-1 . 78 25*02 5. 1876*01 -1.751**02 

-1.7456*00 -7.132E-01 -3.8366*01 

6.6256-01 1.639E-02 5.359**00 

2.3216*00 -1.0676-01 1.358**01 

-4.389E-02 -5.6456-03 -7.003*-01 

1.028**01 
1 .045F*01 
- 1 . 0 16F*00 
-2.1 07F *00 
1. 790F-01 

-4.607**01 
3.111 **01 
1.364**02 
4.015**02 
1.330**00 

-4.4836*00 
-2.7566*00 
-1 -461**01 
-3.092F+01 
5.838 *— 0 2 

-1.2236*02 
2.2586*02 
5.620F+ 02 
-6.9146*02 
2.5226*02 

3.1686*01 
-4.3396*01 
-1.0176*02 
1 .6056*02 
-1 .2286*01 

3.5766*02 

-5.7126*02 

7.900E*01 

5.4626*02 

1.6586*03 

1. 095E*02 
-1.819E*02 
-9,035 E*Ol 
2.2446*02 
3.0546*02 

DERIVATIVES OF INVERSE OF Fft*Q. SOUARFO. 

F IRST n*P I VATt V* -5.9646-05 -2.0736 

SECOND DERIVATIVE 

-4.7786-06 -2. 66 1F-06 = 6.0396-06 

NUMBER 1 
-05 

2.105F-06* 

-1.695*-06 • 

-2 .2946-05* 

-9.1226-06 

-2.471E-06 



OFRTVATIVES OF INVERSE OF FP.FQ. SQUARED, 
FIRST D6R1VATIV* 2.746*-05 8.317F 

SFCOND D* R t V AT I VE 

3. 108E-06 2.65 76-06-= -2.781F-06 - 

NUMBFR 2 
-06 

B. '4226- 07* 

1.6766-06 

2.2 8 BF- 06* 

4.213F-06 

1.211 E-06 



OFRTVATIVES OF INVERSE OF FREQ . SQUARED, 
FIRST OER I VATIVF 9.3506-07 6.267* 

SFCOND DERIVATIVE 

7.8336-08 7.9216-09= -9.4686-08 

NUMBER 3 
-08 

6.346*-09* 

3.0856-08 

5 . 4216-09* 

1.422F-07 

8 . B466-09 



DERIVATIVES OF INVERSE OF FR60. SQUARED, 
FIRST 06RIVATIV* -1.634*-07 -2.658F 

S*COND DERIVATIVE 

-8.6696-09 4.52 86-09 * i.654*-08 

NUMBER 4 
-09 

2.6976-10* 

-3.9266-10 

4.2146-09* 

-2. 4826-08 

4.5166-11 



DERIVATIVES OF INVERSE OF FREQ. SQUARFO, 
FIRST DERIVATIVE -9.6686-08 -2.3B36 

$ *COND DERIVATIVE 

-1 .62 86-08 -4.2556-09= 9.791F-09 

NUM8FR 5 
-08 

2.413F-09* 

-1.1326-08 

-3 .4626-09* 

-1.475F-08 

-3.2056-09 




to 

cO 
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CO 

PROJECTED CROSSING 
l.0310E*01 
1 •3577E*01 
2.7640E*01 
1 • 0 0005*03 
1. 0000E*03 

ITERATION 1 

19. 7500 10.3100 13.5771 27.6*03 1000.0000 1000.0000 


FOUND 2 ROOTS, RFJ FOR THE NEXT ROOT PREDICTED = 27.6403 ,15 BEYOND THF RANGE 


ITERATION NO. 2 DENSITY » 7.8Q61E-08 SQRTI SEA LEVFL DFNSI TY/DENS! TY) * 1.2121E*00 

FLUTTER VEL « 7.6509E*03 AIR SPEED * ?.6504F*03 I VEL-A IR SPEED! *100/ VEL = .0065 


O SUBROUTINE LEGROOT 

DAMPING 
9.3771E-02 
2. 4933 e -0Z 
-B.3046E-02 
-1.9804E-01 
-4. 0967 C -Q1 


FIRST OERIV S c CPN0 ORRIV 

2.2666E-02 4.4264F-03 

3. 731 7F-03 -9.5807E-05 

-7.7608E-O2 

-2.38B7E-02 -4.7313E-03 

“1. 5965 c - 01 -1.0542E-01 


RADICAL IN L c GUERR? EIGENVALUE NO. 
9. 86 T 1 E-"3 5 2 

1.6314E-05 3 

1.1078E-04 4 

-3.6638E-D4 1 

-1.7634E-D2 5 
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APPENDIX B 


FORTRAN PROGRAM LISTING 


The FORTRAN program listing for program MATCH and related subroutines are 
presented in this appendix. 


Program MATCH 


0VFRLAY(MATCH,0,01 

PRD GRAM MATCH { INPJ T-l » OlJTP IJT = 1 , T A , T A 7 = 1 1 T APF8R , 


1 

C * ******* ********** 


c* 

KUMAR 

G- SHAT I 

c* 

FINDS 

FLUTTFR 

c* 

P C RF = 

MAX I MU 

c* 


TOLFPA 

c* 

MAXMAT 

= MAXI 

c* 

I TROPO 

0 p F I NF 

c* 


= 0 , I 

c* 


= ?T I 

c* 



c* 


= -•• » I 

c* 

IMATCH 

= 0 CO 

c* 


OP 

c* 


* 1 MA 


TAP- 5=TNPUT, TAPF6=OUTPU T ) 

**$************** ********* * *** ***** ***** ** ****** ***** 
A, JJLY 24,19 72 * 

MATCH POINT ns CROSSINGS FOR SPECIFIED D P N$ITI C S * 
M (V c lOCI TY-SPFFO of SnUNn*MACH)*lOO/VCLnCTTY * 

NC C SPECIFIED FOR T PR M I NAT I ON (IN PpRC c NTJ * 

MUM NUMBER OF IT r R&TIONS ALLOWED * 

S TH- INHTAI DENSITY FOR MATCH POINT SEARCH * 

NIT I AL DENSITY = S C A LE V C L DENStTY * 

NITIAL O-MSITY = OFNSITY FOR GEOMFTRIC ALTITUDE OF * 
36,?00 FFFT * 

NITIAL D c NS IT Y = «Hn(l) IN NAMELIST F*~>R L p FCROS * 

MPUTES FLUTTfr CROSSINGS AND VELIKI TIPS TOR SINGLF * 
MULT T PL e DENSITIES AS SPFCIFIFD IN l. F FC RQS * 

TCH POINT IS COMPUTED WITH INITIAL OfNSfTy SP p C.IFT c D* 


c* 


BY ITROPO 

c* 

R C F $LD = 

FFFFRFNCF SFA le v c l 0 

c* 

UNI TL = 

1 IF ALL LENGTH UNITS 

c* 


NO. OF LENGTH UNITS / 

c * 



c * 




* 

c 14-SI TY TN APPPOPP.TATC MASS UNITS * 

ARE IN FFFT, NO INPUT R c OU!P c D * 

FOOT, MUST BF INPUT FOR MATCH-POINT* 

search wh c n the length units * 

5 C|_pfjFp arc: OTHFR THAN FF C T * 


Q **************** ************************* ****************************** 


RFAL MACH 

NAMFL I ST /NA MATCH/ P c R F f M A X MAT , MACH , I TRQPP, I MATCH, R FF.S L 0 , UN I T L 
UN! TL = 1,0 
READ 15 ,NAMATCH) 

WRITF (6,NAMATCH) 

IF ( IMATCH • EQ. D ) GO TO 100 

CALL CROSMATI MACH, PFRF, MAXMAT, ITROPO » REFSL 0, UN TTU 
GO TO 200 

100 CALL LFFCPOSI IMATCH, RHOM,RFT MI N, SVPL ) 

200 CONTINUE 
FNO 
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Subroutine CROSMAT 


SUBROUTINE CR OSMAfCMACH* PFPF , MAXMAT, I-TROPO , R C F , UN I Tj_ j 
C OMMON/DFP TVS/ DR c »DMU t OVFLtSOPF,SOMU,SDVPl 
COMMON NM,NMAX,NFI G»NV c f. 

P C AL, MACH 

0 F N T R 0 p = 1*8 36826 P 8 2 
I MATCH = 1 
NOF& = 2 
NL = 2 

NL! - NL - 1 
IF ( IT ROP° ) 6,7,8 
6 R HO = -1*0 
GO TO 9 
? phh = PFF 
GO TO 9 

8 R HO = PFF * 2.9639^-01 
ORFI = -RF I F I *D3 F 

SOP F I - 2.0*0RFI*3RFI/RFI - PF!*RFI*$DRF 
R F I = PFI f DP.FI *OFL F 0 .5 *S FI *OFL*DFL 
GO TO 1 
500 RF TURN 

1000 FORMAT (// ,43H MATCH-POINT IT E D AT ION DIO NOT CONVFRG c TN ,13, 

1 11H l TFR AT IONS ) 

2000 FO°MAT ( / , 7 5H ARGUMENT OF RADICAL IN LAGUFRsr = ,F!2.3,18H, ITFRAH 
ION NO. = ,13, 12H , DENSITY = , C 12.4,8H, V C L = , F9 . 3, / , 10X , 24H , SPff 
20 OF S CU ND*MACH = ,F9.3) 

7500 FORMAT | 1 H 3. | 

3000 FORMAT I//, 14H ITERATION N0*,T3,11H DENSITY = ,F12.4, 35H SORT(SFA 
1 LEY FL DFNSITY/DFNS IT Y 1 = , F12.4,/,16H FLUTTER VFl = ,F12.4, 

213H AIR SPEED = , F 13 *4 , 76H ( V c l-AJ PS PEFD ) * 100/ V C L = ,F8.4) 

C ND 


Subroutine DERVDEN 


SUBROUTINE OFRVOFN {RHHtNDER, RFT } 

C* R HO = SQRTl REFERENCE DENS 1 Ty /OEN S ITY J 

C QMMON/DFRI VS / DRF , 0M!.J,0VF l, SDPF ,SDMU, S OVFL 
COMMON NM,NMAX,N5IG,NVPC 
C0«MC!N/BLK1/ AF,D4F t $0AF 
COMMON /BLK2/ EIG,VEC,AV C C,D!FS,5D!FS 
C0MM0N/BLK3/ TP,TP2 

COMPLEX VSAU, R ! 1 , >U4 ,R 1 5 , R IT ,R I7,RI8 

COMPLEX API 12 , 12 ) ,DAF( 12,121 ,SOAF{12,12) ,E IG U 2 ) , VFC ( 12 , 12 1 , 

1 AVEC ( 12,12),DIFS<12),SDTFS(12},TPU2,12},TP2( 12,121 

C ******** *************************************** 4* ********* ************* 


C* KUMAR G. BHATIA, JULY 2l,!072. * 

C* COMPUTES DERIVATIVES WITH RE$P r T.T TO SORT ( 2 . 378 E-03/DFNS ITY 1 * 

C* D VFL * FIRST DERIV OF VELOCITY, SDV P L = SECOND DFRIV HE VELOCITY * 

C* DMU , SOMU ARE THE FIRST AND S c COND DER!V,RPSP, OF FLUTTER FR* : 0<'*2 * 

C* DRF t SDRF ARE THE FIRST AND SECOND DFR.!V,RFSP, OF REDUCED FR c O. * 

C* NDER * NUMBER OF DERIVATIVES R FOU I p FD» 1 OR 2 * 


C ****************** *******+******^*«** + *«C********<t****!(!«**«t*^3jl ********** 
REWIND 4 

READ {41 M,RFI , VEL , VS AU , E I G , V EC , AVEC , AF , T P 
RHOS * RHO+RHO 
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C IGM = FIG(M> 

CALL T MM PROD ( AVEC, AF , VEC , NM, NV r C , NMA X, NOER »TP2 I 
C IGMI = 1.0/EIGM 
RU = 2. 0*T P? (M, M I /P HO 
R2 = -REAL! TP ( m,M> > 

A2 = - AI MAG I T P( M ,M ) I 
ORF * -A|MAG< RI1 >/ A2 
DMU = EIGMI*( RFALIRI1 1+DRF*R 2 ) 

CQPF = 0.5*EIGM*DMU - PF I * OR F 
DVFL = VEL* CO EF 
PRINT lOGOf RF1, OR F, OMUf DV C L 
IF (NDFR .EQ. 1) RFTURN 
R14 = 2. 0*F I6M*OMU*OMU 
R 15 * 0.0 
9 CONTI NUF 

CALL L FFCROSI I MAT CH, R HO, R F I, VFl I 

OFN = SORT! RFF/RHT I 

O e NL = DFN 

DENI = 1.0/DEN 

GO Tq 2 

1 CONTINUE 

IMATCH = IMATCH+1 

V C L L = VEl 

OFNL = DFN 

RHO = P.FF«0FN!*DFNI 

CALL FOMATCHI I MAT CH, RHO, RF I , V c U 

2 CONTINUE 

CALL SOUND ( MACH, D C N, SOS, D$HS, SDSriS) 

SOS = SOS*UNI TL 

osns = osns*uNiTL 

SOSOS - SOSOS *UN I T L 
F = VEL-SOS 
PPR = F#100 .0 /VFL 
PRINT 2500 

PRINT 3000, I MATCH , RHO,OPN,VFL ,SOS,PFR 
WRITE (7,3000) IMATCH, P.HO, OEN,VFL,SOS,P P R 
IF l ABS(PFR) - PFRF 1 500,500,5 
5 IF (IMATCH ,N E . MAXMAT) GO TO 10 
PRINT 1000, IMATCH 
W RT T E I 7, 1000 I IMATCH 
GO TO 500 

10 CALL OERVOFNI DFN,NDER,RFI ) 

IF ( OFN . NF * OENTROP) GO TO 200 
IF {F .GT. 0.0) GO TO 200 
D SO S = 0.0 
SOSOS = 0.0 
200 CONTINUF 

OF * DVEL - D SOS 

SDF s SDVEL-S DSOS 

H = NL1*(NL1*DF*DF-NL*F*SDF) 

IF (H ,GE. 0.0) GO TO 250 
PRINT 200C, H, IMATCH, RHO, VELtSOS 
WRITE ( 7,2000) H, I M ATCH , RHO , V El , SOS 
GO TO 500 
250 CONTINUE 

H = SQRTCH) 

H * S I GN (H, DF ) 

OEN * DFN - NL*F/(DF«-H) 

IF IDENL - DFNTROP) 260,290,270 
260 IF (OEN - DFNTROP) 290,290,280 
2^0 IF COEN - DFNTROP) 280,290,290 
280 OEN = DFNTROP 
290 CONTINUE 

DENI = 1.0/nFN 
OEL * DFN - D EM 
00 100 L = 1 » NV EC 

IF (L .EC. Ml SO TO 100 

RIT = DRF*DRF*TP(L ,MI *tp(m,L ) - 2.0/ RHO *DRF*< TP2( L » M) *T-p ( M,L ) ♦ 
1 TP2(M,L )*TP(L,M)I «■ A . O/RHOS+T P2< L , M ) *TP 2 ( M, L ) 

RIT = RIT /( 1.0-E IGM/E IG( L) > 



APPENDIX B - Continued 


«I5 = r i 5 + R i t 
100 CONTINUE 

R 15 = 2.0*EIGMI*RI5 

RI7 = -E IGM I * ( 6. 0*TP?( M, /»H0S-4.O*TP (M,M1 *DRF/RHn*VS AU*0RF* 

1 DRF I 

RI8 * -EIGMI*TP(M f M) 

SD&F * -AIMAG|RI5*«I4+RI71 /ATMAGIRI8) 

SOMU = RE AL ( P 15* Rl 4+RT 7) ♦ R EA L ( R 1 8 ) +SDRF 

SDVEL = DVEL*COFF + VEL*( 0»5* P IGM«(SDMJ- C I GM«OMU*DMU I * PFTMRFI* 
1 DRF*DRF-SDRF) ) 

P RI NT 2000, SPRF, SOMU.SDVFl 
RETURN 

1000 F HRM AT ( / / ,1 7H REDUCED FRFQ . = ,Ffl.3,/,63H **THE DFRIVATIVES AR C W. 
1R.T. SORT (S FA LEVEL DENS IT Y/ DENS ITYl **,/,/, 28H DEPIV. OF REDUCED 
2FRFQ. * t F 13» 4, / , 2 2H DERIV. OF FREg**2 = ,E13.4,/, 

3 2 3 H D C P I V* OF VELOCITY = ,^13. 4) 

2000 FORMAT!/, 22H SECOND OFRIV. OF RF =,E13.4,29H , SECOND D C PIV. OF FR 
1 C Q**2 =, F13.4,24H , SECOND OFRIV. OF VFL=, C 13.4) 

END 


Subroutine SOUND 


SUB POUT INF SOUNOIM ACH, PHO, SOS* OSOS t SOSO S ) 
REAL MACH 

C* RH3 = SQRT(S p A LEVEL OENSI TY ) / SQP T ( D*= NS ! T Y ) 
IF (RHD .GT. 1 • 836 826 B82 ) GO TP 10 
A = 1515.639571 
3 = -520.6920622 
C = 121. 182 A9 16 
SOS = MACH* ( A+RH0*( B-*- p HD*C I ) 

OSOS = MACH*I B+2 • 0* R H0*C ) 

SOSOS = MACH * 2.0 * C 
R F TU °N 

10 SOS = MACH*96B.08 
OSOS = 0.0 
SOSOS = 0.0 
R c T U R N 
c NO 


Subroutine LEGROOT 

SUBROUTINE LFGPPOT { N c T G , I C ON , P F I , G , I AR , I PRT , NL ) 

C0MM0N/RLK2/ c IG,VFC, AVFC, 01 FS, SDTFS 

COMPLFX P I G 1 1 2 ) » V F C ( 12,12) , A V e f ( 12*121 , D I F $ ( 12) , SDIFSU. 21 

0 I^F NS I ON PFI (1) ,1 CON lli ,G(lMAPm 
REAL I OK » M 1 NI OK 

C* KUMAR G. BHATIA t JUNE 12, 1972. * 

C* COMPUTES THF ROOTS USING NOTIFIED LE GU C RPE ITERATION, WHERE THE * 

C* ROOTS CORRESPOND TO THE IMAGINARY PART OF INVFRSE OF THF FRFQ.* 

C* SQUARFD AS A FUNCTION OF INVERSE OF REDUCED FREQ.(RF). * 

C* AT INPUT F: F 1 1 11 CONTAINS 1 /RF WHEOC PUNCH ON+DER I VATI V c $ ARE KNOWN 

C* AT OUTPUT PFI I J) CONTAIN PROJECTED ROOTS * 

C* PROJECTED CROSSING = 3000, R FA 1_ PART OF F I GENV ALIJE IS N^GATlv* * 

C* PROJECTED CROSSING = 2000, 1C0NIJ) . NE . 0 * 

C* PROJECTED CROSSING - 1000, REAL RFKJJ COULD NOT B c PREDICTED * 

C« **$**##«*#*** &**** 4 4 4 4: ****** **** ******* ***** ** ******************** 

1 OK = RFItl ) 

RFI2 = 0. 5*1 OK 
N = 0 

CALL DAMPAR (NEIG, E IG.G, IAR » 
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IF I IPRT .NF. 01 PRINT 1000 
DO 200 J -1 * NF IG 
I = IARlJl 


II = I 

G l = 0.0 
G 2 = 0.0 
A ~ 0.0 

IF I GUI *F0. -1000.0 I GO TO 4 
IF IICCNIJ) * EQ . 0) GO TO 5 
RFIIJI = 2000 .0 
RFIJ = RFIIJI 
GO TO 100 
RFIIJI = 3000 .0 
RFIJ = RFIIJI 
GO TO 100 


5 C n NT I NUF 

&0 = AIM AG ( FIG ( I I I 
A 1 = AIMAGIOIFSini 

A3 = AIMAGISDIFStl) 1 >( .„ fCC1 

RO=RFAUEIGUU » Rl-REALlOIF5(III * P2=RFAU I SOIFSI 
GO * GII) $ G1 = (Al - GO*RU/RO S G2 = {A2-G0*R2- 

A a G1 *G 1 - G 0*G 2 _ 

IF (NL .NF. 0) A = (NL-l)M (NL-1 I*G1*G1-NL*G0*G2 ) 

IF ( A ,GT. 0.0 I GO TO 10 
IF (GO .LT, 0.0) GO TO 8 
IF ( G 1 .GF. 0.01 RFIIJI = I OK - G0/G1 
IF ( G1 .LT. O.Ol RFIIJI = IOK + G0/G1 
IF (RFIIJI .LT. RFI2) RFIIJI = RFI2 
GO TO 12 
8 CONTINUE 


N = N + 1 

RFI IJ I a 1000.0 
RFIJ = RFIIJI 
GO TO 100 

10 IF (NL .EQ. 0 > GO TO 11 

RFIIJI * IOK - NL* GO/IGl+SQRTI A) J 
GO TO 12 

11 RFIIJI = IOK - GO/ SORT I A I 

12 continue 

RFIJ = RFIIJI 

IF I J .EC. 1 I GO TO 100 

IF I RFIJ . GF , RFIIJ-1J I GO TO 100 

Jl * J - 1 

on 15 1*1, Jl 

IF ( RFIJ . GF . RFIIJ-II ) GO TD 16 
JMIN = J - I 


15 CONTINUE 

16 SWF! = IAR ( J ) 

J JMIN = J“" J MI N 
DO 20 I=1,JJMIN 

I AR C J-I *1 > = I AR C J -I I 
20 RFI ( J-I+l) = RFI U-I) 
RFI ( JMIN) = RFU 
IAR l JMIN) - SAVE I 


100 CONTINUE 

IF { I PRT ,F0. 0> 30 TG 200 
PRINT 1500, GO »G1 ,G2,A, I I»RFI J 


200 CONTINUE 

IF I RFim .GT. 0.0 ) GO TO 400 
IF ( RFI(NEIG) .LF. 0.0 ) GO TO 400 
NFIG1 * NF I G - 1 
00 300 J = 1, NF IG1 

JJ - N c I G1 - J F 1 
IF (RFKJJM 250,250,300 
JJ1 * JJ * 1 
DO 300 I =1 ,JJ 


1 1 ) 

2.*G1*R1)/R0 


250 
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RF I U) * RFI ( JJ1 ) 

300 CONTINUE 
400 CONTINUE 

IF I N . LT. N F I G J RETURN 
PR. I NT 2000 , N 

1000 FORMAT {/ t 19H SU8RD UT INF L F GROOT, / , 6X , 8H DAMP TNG, 10X , 1 2H FIP$T OFRI 
1V,7X,13H SFCCND DERIV,4X,20H RADICAL IN LFGUFR RF , 3X , 1 5H r I G C NV At U p 
2 N0..,3X,19H PROJECTED CROSSING! 

1500 FORMAT (4 ( 5X,E 11.4, 5X ), 10X, I 2, 3 3X, C 11 .4) 

2000 FORMAT (//,38H ARGUMENT OF LAGU c RR c IS NFGATIVF FOR ,I3, 7 h MODES ) 
C ND 

Subroutine DAMPAR 


SUBROUTINE D AMPAR ( N c 10, E I G, G» IAR) 

COMPLEX FIG tl ) 

DIMFNS ION Gm.IANlI 
[ARdl = 1 

GO.) = AIMAGI fig(I) J/RFAL < C IG( 1 ) ) 

IF IRFAL( C IG( 1M .L c . 0.0) G(l) = -1000.0 
on 5 I -2 *NF IG 
I AR II) = I 

G ( I ) = AIMAG(FIG( I) )/R fr AL( c IG(I )) 

IF (REAL (F IG( I) ) .L p . 0.0) GUI = -1000.0 
IC - I 
11 = 1-1 
DO 4 J = l,I 1 

M = I 1 -( J-l ) 

ICC = IAR( IC) 

MM -= I AR ( M ) 

IF (GUCCI . L c • G(MM|) GO TO 5 
IT = I A R ( M ) 

IAR(M) = IAR 1 1 Cl 
IAR(IC) = IT 
IC = IC-1 

4 CONTINUE 

5 CONTINUE 
RETURN 
END 


Subroutine GETAERO 

SUBROUTINE GETAERO (NM.NMAX ,RFI , IDiRF IL,REIR t OEL> 

C0MM0N/BLK1/ AF, DAF,SDAF 

COMPLEX AF( 12,12) t 0AF( 12 , 1 2 ) , SOA F ( 12 , 12 ) 

C* GETS AERODYNAMIC FORC c S FROM R AN DOM ACCESS FIL C 88 IF ID=0,FLSE * 


C* GETS DERIVATIVES TOO. ENTRY GETDAER G FTS DEPIVS. ONLY. * 

C* REI = 1.0/REDUCED FREQ., RFIMIN = FIRST PFI RECORD ON 88. * 

C* DEL = CONSTANT INC REM r NT OF RF I ON 88. * 

C ^ KUMAR G.9HATIA, JUNE 13,1972. * 

C* 

IF ( RF I .GE. RF I L .AND. R F I .LF. R F I R ) GO TO 10 
PRINT 2000, RF I 


2000 FORMAT!/,* RF I = *,F10.3,*, IS OUTSID c THF RANGF OF VALUES *) 
STOP 
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10 CONTINUE 

NW = 2*NMAX*NMAX 
STEPS = ( RF I- RFI L) /DEL 
IK = STEPS 
I R - 2 

IF ( ( STEPS- I K I • LE . 0. 5 > IR=1 
PFI - RF IL ♦ 1 1 K ♦ I P-1)*DFL 
IK = < IK4- IR-1 ) *3 ♦ 1 
CALL READMS(88,AF,NW,IK) 

IF (ID .EG. 0) RETURN 
ENTRY GETDAER 

CALL READMS(88»DAF,NW t IK«-l) 
CALL READ-M${88,SDAF,NW, IK+2I 
R C TURN 
F NO 


Subroutine RAND AX 


C* 

C* 

C* 

C* 

C 

c 

c 

c 

c 


S UR POUT INF RANDAXI NM ( NMA X, SS » RR i PHO, NRF ) 
dimension NRF U I 
C0MM0N/BLK1/AF ,D AF t S OAF 

COMPLEX AF( 12 « 12) • OAF ( 12, 1 2 ) » SOAFt 12 , 121 
DIMENSION R0(2»12»121»Rl(2tl2*12ltR2*2tl2»12) 
f OU I VALENCE ( AF , RO ) * ( OAF t R 1 ) » t SDAF , R 2 ) 

Thf SUBROUTINE READS FROM TAPF A AND* TP ANS FER S TO RANDOM ACCESS * 
FILE ON TAPE 88. THF AFRO FORC.FtOERlV MATRICES ARE MULTIPLIED* 
BY DENSITY PARAMETER B ?F QR F TRANSFFR TO g8. * 

KUMAR G.BHATIA, JJNF 13,1972 


COMPUTE THE DENSITY PARAMETER 
DP IS IN LB.SEC**2/INCH UNITS 

INPUT SS= i * * * S SEMlSPAN»BR=REFE PENCE SFMICHORO ARE IN INCH C S 
RHO= AIR DENSITY IN SLUGS/FT**3 
PI * 3.14159265358979 

DP = 4.0*PI*BR*SS*SS*RH0 
REWIND 4 

R EA D (4 ) NK»MACH,NM 

CALL 0PENMS(88»NRF ,1600, 0) 

NW = 2*NMAX*NMAX 

DO 100 I K-l * NK , 4 , , 

READ (4 ) RF, X , ( ( R 0 < l,I,J)fI=l,NM) »J=1 , NM t » C (R0(2,!»J),I=lfNM) , J-l f 

1 NM> 

RF A D 14 ) ( IR1 ( 1 , 1 1 J ) f 1 = 1 1 NM ) , J= 1 , NM), ( (R1 I 2 » I » J ) » I - 1»NM ) , J=1 * NM ) 

R C AD {4 ) U R2 ( 1» I f J 1 » 1= If NM > t J*1 ,NM 1 , MR 2 ( 2 » I » J ) * I - 1»NM) » J = 1 » N M ) 
RFI = 1.0/RF 
F = DP*RFI«RFI 
DO 10 1=1, NM 
DO 10 J*l, NM 
AF ( I, J) = F*AF { I, J) 

OAFIIfJ) = F*DAF(I,J> - 2.0*RFT*AFU, J) 

S DAF I I , J ) = F*S OAFl I, J ) - 2 .0 *RFI *( 2 . 0*DAF ( I , J ) + RF T *AF( I , J )) 

10 CONTINUE 

IK3 = I 1 K— 11*3 * 1 

CALL WPITMS188, AF,NW, IK3) 

CALL W D ITMSI88, DAF,NW, TK3 + U 
CALL WRITMS(88, SDAF,NW, IK3 + 2) 
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100 CONTINUE 
RETURN 
END 


Subroutine LEFCROS 

S UBRDUT I N c LEFCROS (TNftTCH, RHOvi t REIMIN,SV c L » 

COMPLEX AF(l2tl2l»D4F(12,12)»SnAFC12fl2l,TPI12tl2),TP?«l? f 12], 

1 EIGC1 2),VEC(12 f 12) *AVFC(12,I 21 rOIFSH2» , SOIFS'd 2 I 

DIM C NS ION I NTH (12, 2) , NR F (i 600 ) , I CPN( 12 ) , SM ( I 2 , 1 2 I , SK< 12 , 12 ) , 

1 C ( 12,12) ,RFI(I2J ,G( 12) ,IAR (12) ,VEU 12) ,RHO( 10) 

EQUIVALENCE (SK,C) 

C DM PL C X MUM 

COMMON NM,NMAX,NEI G,NVEC 
COMMON/BLKl/ AF, OAF, SDAF 
C0MM0N/BLK2/ FIG , VEC , A VEC , 01 FStSDTFS 
C0MM0N/BLK3/ TP,TP2 
C0MM0N/BLK4/ 5M,C, INTH 

NAMELIST/NAM1/ SK , SM, L STI F F , $S , B P , NM , RF 1 1. , RF I R , DFL , 

1 NRODT »N!TMAX,NP»RHO,RFIMI N, I PRT , IOPT 

NAMFL I ST/OPTION/ N MAX , NE IG, N EVRF D, FI GRAT,NL, ICON 
C* *********4**************** ** ##** ****** 4* *************** ** ***« 

C* COMPUTES FLUTTER CROSSINGS AND VELOCITIES FOR SINGLE OR MULTIPLE* 
C* DENSITIES, IMATCH * 0 AND NO = NO* OF D C NS!TIF$. * 

C* FOR IMATCH *NE. 0, THE LOWEST FLUTTER VELOCITY AND OTH c R INFO IS* 


C* RETURNED TO THE CALLING PROGRAM. THE INITIAL GU C S$ F0p PFIMIN* 

C* IS PICKED UP FROM THE NAMELIST FOR IMATCH = OOP 1, FOP OTHER* 

C* VALUES OF IMATCH THE GU^SS SUPPLIED FROM TH C PAPAMFTFR LIST. * 

C* KUMAR G. BHATIA, PROGRAM CHECK COMPL ET C D JULY 20 , 1972 . * 

C* ****************«*+««*««««^**J(:** *********««♦♦#♦♦♦**♦;**♦* fr***^:**** 

C* DEFINITION AND ASSIGNEM c NT OF THF NAMFL I ST NAM1 PARAM r T FPS * 

C* LSTIFF * 0 DIAGONAL STIFFNESS MATRIX IS INPUT IN SK, DIAGONAL * 

c* flexibility matrix is computed and STORED !N SK * 

C* LSTIFF = *1 FULL STIFFNESS MATRIX IS INPUT IN SK, IS INVERTED AND* 

C* DESTROYED USING COC matrix INVERSION POUTTN c * 

C* LSTIFF * -1 DIAGONAL OR FULL FLEXIBILITY MATRIX IS INPUT IN SK * 

C* SK = GENERALISED STIFFNESS OR FLEXIBILITY MATRIX, SE C LSTIFF * 

c* sm = generalised mass matrix, maybf diagonal or full * 

C * SS * SEMISPAN, RR = REFERENCE CHORD - BOTH MUST BE IN APPROPRIATE* 

C* UNITS * 

C* DEL = EQUAL INCREMENT I1N PFI AT WHICH AF PDDYNA MI C FORC r $ ARE ON *WP 

C* TAPE A * 

C* RFIL = MINIMUM VAUJP OF RFI F n R WHICH AERO FORCES AR C SUPPLIED * 

C* R.FTR ■= MAXIMUM VALUE OF RFI FDR WHICH AFRO FORC c S AR C SUPPLIED * 

C* TAEPO * 0 AERODYNAMIC FORC c S AND FIRST TWO CFBIVATIVFS APF * 

C* SUP PL I FD AT EQUAL RFI INTERVAL OF DEL, STARTING WITH RF U 

C* IAFRO .NE. 0 ONLY AERODYNAMIC FORC c S ARE SUPPL^O FOR INCREASING * 

C* VALUFS DF RFI, STARTING WITH RF T L * 

C* NM = NUMBER OF MODES, NM AX = MAXIMUM NO. OF MOD c S ALLOWED =12 * 

C* NEIG = NO. OF EIGENVALUES TO BE COMPUTED, N'^IG .LE. NM * 

C* NV e C - NO. OF EIGENVECTORS TO R c COMPUTED, NVEC .LE. NM * 

C* FIGRAT REQUIRED ONLY WHEN NEVRED .NE. o, SFF NEVRED * 

C* NEVRED = 0 NEIG EIGENVALUES AND NVFC VECTORS ARE COMPUTED FOR * 


C* FIRST AND SUBSEQUENT C I C,fn SOLUT IONS * 
C* NEVRED .NF. 0 AFTER TH C FIRST E IGENSOLUT I ON ONLY THE SMALLEST * 
C* EIGENVALUES (AND VECTORS) ARE COMPUTED SUCH THAT THE SMALL P ST * 
C* c IGENVALUE NOT COMPUTED IS AT LEAST FIGRAT TIMES T H c FLUTTER * 
C* EIGENVALUE * 
C* ICON IS INITIALLY SET TO ZERO, IF ICQN(L) IS INPUT AS NONZFRD * 
C* L TH LARGEST DAMPING ROOT PROJECTION FnR F I NO TNG PFI, IS * 
C* NOT CQMPUT ED * 
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C* NITMAX = MAXIMUM NUMR C R OF ITERATIONS ALLOWED PER ROOT 
C* NROPT = NO. OF ROOTS TO BE SEARCHED * 

C* NO = NO. OF DENSITIES FOR WHICH THE FLUTTER CROSSINGS A R c DE S I R C D* 

C* RHO IS THF VECTOR OF DENSITIES IN APPROPRIATE MASS UNITS * 

C * IPRT = 0 SUMMARY PRINTOUT ONLY, -l ROOT PROJECTIONS AMD EIGFNVALU* 

C* E S PRINTED AT F AC H STEP, * 2 EIGENVALUE D C RIVS. AND EIGEN* 

C * VECTORS ALSO PRINTED AT EACH STEP 

C* NL SPECIFIES THE ASSUMED NO. PE ZEROS OF DAMPING AS A FUNCTION OF* 
C* RFI, NL=0 ASSUMES DAMPING AS A T R ANS CE DENT AL FUNCTION, DFFAULT * 

C* V AL UF IS 0 * 

C* RFI = RECIPROCAL OF R c DUCF D FP EQU C NC Y * REDUCED VELOCITY * 

C* RFIMIN = INITIAL GUESS FOR RFI * 


REWIND T 
DO 5 1-1,12 
5 ICON! I ) =0.0 
NL = 0 
I OP T = 0 
IPftT * 0 


N? V P C D = 0 
NMAX = 12 
PEADI5,NAM1) 


WRITE! 6,NAM1) 

N c I G - NM 
NVEC = NM 

I e ( TORT ,EQ. 0) GO TO 7 
READ 15 , OPTION 1 
WRITEI6, OPT ION) 

7 CONTINUE 

DEL 2 = DFL/2. 0 
IF ILSTIFF) 14,6,12 
6 00 10 1=1, NM 
DO 10 J s 1 » N M 

IF I J .FQ. I ) C(ItJ)=l*0/SK(!,J) 

1.0 CONTINUE 

GO TP 14 t _ 

12 CALL MATI NV ( SK ,NM , DUMMY, 0, O c T,NRF,INTM,NMAX,ISCALE) 

14 REF RHO = RHOt 1) 

IF IRHOM .EG. -l.Jt RHOM = RH0I1) 

IF < IMATCH .NE. 0» REF RHO = PHPM 

15 CALL RANDAX(NM,NMAX,SS,BP,RFFRHO,NRFI 
ENTRY FOMATCH 

IF (IMATCH .EQ. 0) GO TO 16 
ID = 1 

RHn ( ID) = RHOM 

16 CONTINUE 

DD 500 TD=1 ,ND 
PRINT 9000 

PRINT 4000 , RHC1U0) , RFIMIN 
W RI T E ( 7, 40D0 ) RHO( ID),RFIM IN 
NR = 1 
20 CONTINUE 

DO 100 1=1, NITMAX 
PRINT 9000 

CALL GETAERO( NM, NM AX ,RFI MI N,0,RFIL,RFIR,D c O 
IF (ID .EQ. X .AND. IMATCH . LE . 1) GO TO 40 
DM = RHOUDI/REFRHO 
DO 30 I A = 1 » NM 
DO 30 JA = 1» NM 

30 AF 1 1 A, J A ) = DM * AF(IA,JA) 

40 CALL EIGSDLINM, NMAX, NFIG, NVEC, IPRT) 

CALL GET DAERINM, NMAX, RFIMI N, l , RF I L ,RF I P , DEL ) 

IF (ID .FQ. 1 .AND. ! MATCH .LF. 1) GO TO 60 
DO 50 I A = 1 » NM 


DO 50 J A=1 , NM 

OAF ( IA, J A ) = DM * DAF ( I A , JA I 
50 SDAF ( I A , JA ) = OM * SDAF(IA,JA) 

60 CALL DERF(2 »NM, NMAX, NVEC, RFIMIN, IPRT ) 
R.FH1) = RFIMIN 

CALL LEGROOTl NFIG, ICON,RFI ,G, 1AR,IPRT,NL) 


39 



APPENDIX B - Continued 


PRINT 1000, I. RFIMIN, <R F It J > t J* 1 , N c IG 1 
IF INFVRFD .FQ. 01 GO TO 80 
NFVRFD * 0 
I M » IAR(NR) 

MUM * FIG(IM) 

I Ml « TM ♦ 1 
DO 70 J*IM1,NM 
0 * MUM/F1G(J> 

IF (D .GE. FIGRATI GO TO 75 
70 CONTINUE 
NEIG * NM 
NVEC « NM 
GO TO 80 
75 NET G * J-l 
NV r C « J-l 
PRINT 6000, NVEC 
WRIT F( 7, 6000) NVE: 

80 IF ( ABS(RFIMIN-RFKNR)) ,L C . DEL2 ) GO TO 2 10 
RFIMIN * RFI(NR) 

IF (RFIMIN .GT. RF IL .AND. RFIMIN .LT. RF I R ) GO TO 100 

NR1 * NR - 1 

PRINT 7000, NR1, RFIMIN 

WRITEI7, 70001 N«l, RFIMIN 

IF ( NR1 .GT. 0) GO TO 400 

STOP 

100 CONTINUE 

PRINT 1500, N R, I 
STOP 

110 CONTINUE 

IM « IAR(NR) 

PRINT 5000, IM, (EIGII), 1=1 , NEIG) 

PRINT 5500, U AR( I I , I«1 ,NFIG) 

VEL(NR) = BP*SQRT( 1.0/RFAL (EIG(IM)| J*RF!MIN 
IF I NR .NF. 1 > GO TO 200 
120 CONTINUE 
R C W I ND 4 

WRITE! 41 IM ,R FIMIN ,VFL( NR I , T P2 (I M , I M I , F I G , VFC , A VFC , AF , TP 
NFRQOT * NR 
$V P L « VFL(NR) 

200 IF (VELINR) .GE. V EL INFROQT ) ) GO TO 300 
GO TO 120 

300 PRINT 3000, N R , VEL ( NR) , RF I M! N, I 
WRITE (7,3000 I NR , V EL (NR I , RF I M I N, 1 
310 CONT INUF 

IF l NR . FQ. NROOT ) GO TO 400 
NR * Nfifl 

IF ( ABS(RFIM IN-RFKNR >1 .GT. DF12 > GO tq 350 

RFI(NR) * RFIMIN 

WRI TF ( 7, 80001 RF I ( NR ) » DEL2 

PRINT 8000, RFI(NRI,DFL2 

TM « I AP (NR > 

VEL ( NR I = 8R* SQRT( 1. O/RFAL ( c |G ( I M) ) )*RF1MIN 
IF (VELINPJ .GF. V Ft I Nf ROOT ) ) GO TO 320 
REWIND 4 

WRITF( 4) 1M.RFIMIN ,VFL(NP I ,TP?(I M,IM) ,FI G,V P C, A VFC ,AF,TP 

nfroot = nr 

SVFL * VFUNR) 

320 CONT INUF 

PRINT 3000, NRtVELINRI tRFMlN, I 
WPI TFl 7, 3000) NR v VFUNR), RFIMIN, I 
GO TO 310 
350 CONTINUE 

RFIMIN * RFI(NR) 

IF (RFIMIN ,GF. RFIL .AND. RFIMIN ,L F . R F I R ) GO to 20 

NR1 * NR - 1 

PRINT 7000, NR1, RFIMIN 

WRI T P ( 7 , 7000 ) NR1, RFIMIN 
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.FQ. 01 RFIMIN = RFl(l) 
NF • 0) GO TO 600 


500 

600 

1000 

1500 


400 CONTINUE 
IF ( IMATCH 
IF! I MATCH , 

CONTINUE 
CONTINUE 
RETURN 

FORMAT (// 10 H ITERATION, 13 , /,F 1 ?. 4 , 8 X, 5 F 12 . 4 ,(/ 20 X, 5 F 12 . 4 ) ) 
FORMAT!/,* PROGRAM TEPMINATFO, COULD NOT FIND ROOT NO. *,I 3 » * IN* 
1 , 13 ,* ITFRATI ONS* ) 

3000 FORMAT!//, 12 H ROOT NUMBER ,I 3 , 14 H , V c LOCITY = ,F 9 . 3 , 

l 9 H , RF I S , F 8 • 3 , 29 H , NO. OF ITERATIONS R c 0 n, = , 12 ) 

4000 FORMAT!///, 11 H DENSITY * ,^! 5 . 6 , 12 H ♦ PFIMTN = ,F 12 . 4 ) 

5000 FORMAT!/, 26 H FLUTTER EIGENVALUE NO. - ,I 3 , 13 H, EIGENVALUES,/, 

1 ( 1 0 El 2 .4/ ) ) 

6000 FORMAT!/, 91 H OPTION TO REDUCE NO. OF FIG C NVALUES AND c I G P NV C C. T FIRS 
1 FXFRCISF 0 , NFIG AND NVFC SFT EQUAL TO,T 3 ) 

5500 FORMAT !/, 62 H RF I FOR PREDICTED CROSSINGS CORRESPOND TO EIGENVALUES 


1 NU M8 C RS , ! / 12 15 ) ) 

7000 FORMAT!/, 6H FOUND, 13, 42H ROOTS, RF I FOR THF NEXT ROOT PREDICTED = 
1,F10.4,21H ,IS BEYOND THE RANGE) 

8000 FORMAT!/, 35H RFI PREDICTED FOR THE NF XT °OOT = ,F10.4, 

1 52 H, DIFFFRfnCF FROM RFI FOR. PREVIOUS ROOT IS L C SS THAN,F8.4) 
9000 F0RMATC1H1) 
c NO 


Subroutine DERF 


SUBR OUT I NF DE RE { NO, NM, NM AX , NVFC, RFI, IPRT | 

C0MMON/BLK 1 / AF , DA F , SDAF 
COMMQN/0LK2/ EIG.VEC, AVEC, DI FS ,$DIFS 
COMMON/ BLK3/ TP,TP2 

COMPLEX AF( 12 ,12) , DAF! 12 ,121 , St>AF( 12 ,1?) ,EIG(12),VEC(12,12), 

1 AVEC! 12,12 ), DIPS (121 ,S0IFS112),TP(12,12> ,TP2( 12, 121, A,B, C, 

2 MUM 

C ******** * ** * * *** ** * **** ** ** *** ********* ***** * ************* ************* 

C * * 

c* ND = NUMBER OF DERIVATIVES REQUIRED, 1 OP 2. * 

C* NM = NUMBER OF MODES. * 

C* NMAX * MAXIMUM NUMBER OF MODES DEFINING SIZE OF VARIOUS ARRAYS. * 

C* NV C C = NUMBER OF FREQUENCIES FOR WHICH THE DERIVATIVES COMPUTED. * 

C* pp = REDUCED FREQUENCY * 

C* c IG [ NMAX ) VECTOR OF INVERSE OF FREQUENCIES SQUARED. * 

C* VFC ! NMAX , NMAX ) ARRAY OF EIGENVECTORS, ONE PER COLUMN. * 

C* A VEC ( NMAX ,NMA X) ARRAY OF ASSOCIATED EIGENVECTORS, 0N P PER COLUMN.* 

C* AF (NMAX, NMAX l AIR FORCE MATRIX. * 

C* DAF ( NMAX, NM AX ) FIRST DERIVATIVE OF AF W.R.T. R FOUCFD FREQUENCY. * 

C* SOAF (NMAX, NMAX) SECOND DFRIV. OF AF W.R.T. REDUCED FREQ. * 

C* OIFSINMAX) FIRST OERIV. OF *=IG<NMAX) W.R.T. 1/RF. * 

C* SDIFS(NMAX) S,ECONO DERIV. OF EIGINMAX) W.R.T. 1/RF. * 

C* TPINMAX, NMAX) ,TP2( NMAX, NMAX) TEMPORARY STORAGE ARRAYS. * 

C* COMPUTES DERIVATIVES OF INVERSE OF FREQUENCY W.R.T. 1/RE0UC C D * 

C* FREQUENCY WHERE LAMBDA* FR C Q**2 IS DEFINED 9Y THE EQUATION * 

C* (STIFFNESS - L AMBDA ( MA SS *AE ) )VEC =0 * 

C* KUMAR G. BHATIA, JUNE 8,1972. * 

C ********************************************************** ************* 

C* * 

C COMPUTE THE REQUIRED ELEMENTS OF TRIPLE PRODUCT MATRIX 

C AVEC TP A SPQSED * DAF * VFC. 

CALt TMMPROD ( AVEC, DA F , VEC , NM , NV P C , NM AX, 2 ,TPI 
C COMPUTE THE' FIRST D C RIVATIVES OF !./FREQ**2 W.R.T. 1/RF 

RF = 1.0/RFI 
DO TO M= 1 ,N VEC 

10 DIF SIM! * -RF*RF*ETG!M)*TP<M,M> 

IF ( ND * C Q • 1 > RETURN 
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C COM PUT F THE SECOND DERIVATIVES GF l./FREQ*T2 W.R.T. 1/«F 

CALL T MM PROD ( AVFC, SDAF,VEC ,NM,NVFC .NMAX , 1 ,TP2 > 

RF4 * RF**4 
DO 100 M= 1 f NV EC 
MUM = FIGtMi 
A = -2.0*«F*niFS(M> 

B = 0.0 
DO 20 L - 1, NM 

IF I L . EO. M ) GO TO 20 

R = B ♦ TP(L , M I*TP ( M »L I / (1.0 - MUM/F!G(LM 
20 CONTINUE 

B = - 2. 0*B *R F4*MU M 

C = RF4*MUM*TP2 (M,M) 

SnlFS(M) * A ♦ B ♦ C 

IF IIPRT ,CQ. 2 ) PRINT nO'3,H, 01 FS ( M I , S D I F S ( M J ,A,B,C 
IOC CONTINUE 

1000 FORMAT!// 48H DFRIVATIVES OF !NV C RSE GF FR^Q. SQUAREO, NUMBER, 12, 

1 /, 17H FIRST DERIVATIVE, F15. 3, E13. 3, /,18H S c CCND DFR I V AT T V c 

2/E13.3 ,E12.3,lH=, e 13.3t c 12 ,3,1H+,F13 .3, c I2.3,lH+,Fl3.3,c 12 .3) 

RETURN 

END 

Subroutine TMMPROD 


SUBROUTINE TMMPROD (A ,0,V,N ,NV, NMAX ,ND,R) 

COMPLEX A { NMAX *1 I * D ( NMAX , I ) , V ( NM AX, 1 ) , R ( NM AX , 1 ) , TEMP 
C T 

C COMPUTES A *D*V - R . IF N0=1 TH C N ONLY DIAGONALS ARF C.OMPUT FD 

DO 100 1=1, NV 
DO 100 J = 1 » NV 

IF (NO .EQ. 1 .AND. I . NF . J) GO TO 100 
R ( I » J ) * 0. 0 
DO 50 K =1 , N 
TEMP = 0.0 
DO 40 L= 1 , N 

40 TEMP = TFMP ♦ D( K, L > *V l t ,J) 

50 R(I,J> = R(I,JI ♦ A(K,I)*TFMP 
100 CONTINUE 
RETURN 
FNO 


Subroutine EIGSOL 


SUBROUTINE E1GS0U NM, NMAX, N C IG,NVEC, IPPT) 

C0MM0N/BLK1/ AF,H,HL 
C 0MM0N/BLK2 / E IG , V FC , AV FC, CNT, COLM 
C0MM0N/BLK3/ TP,TP2 
C0MM0N/8LK4/ $M, C * INTH 

COMPLEX AF( 12 V 12),H (12,12) ,HL (12,12) ,EIG1 12 I ,V C C ( 12, 12)', 

1 AVECt 12,12 >, CNT [12 ),C0LM(12> ,TP ( 12 , 1 2) , TP2( 12 , 1 2 ) , SUM , SUM1 

DIMENSION SM( 12,12 > ,Ct 12,12) , INTH( 12 ,2) 

T 

COMPUTE THE PRODUCT C*(SM+AF ) 

DO 10 1=1, NM 
DO 10 J= 1 ,NM 

TP( I, J) = 0.0 
DO 5 K= 1 * N M 

5 TP ( I , J) = TP ( I , J ) f C(I ,KI*(SM(K, J >*AF(J,K) ) 
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10 CONTINUE 

INTHU,1) = NM 
I NTH ( 2 , 1 ) = NVEC 

CALL EEC M (T P» E IG, AVFC, HL, H , CNT , C OLM, I NT H » N M A X ) 
IF ( I NTH l It 1 1 .EQ. NM 1 GO TO IF 
PRINT lOOOt INTMltl* 


C 


C 


C 


STOP 

15 CONTINUE 

COMPUTE TP2 * $M*-A F , AND C*<$M+AF) 

DO 20 1 = 1, NM 
DO 20 J=1 » NM 

20 TP2II,JI * SMU, J » +AF< I* J » 

DO 30 1= 1 ,NM 
DO 30 J=1,NM 

TPtltJ) = 0,0 
DO 25 K* 1, NM 

25 TP 1 1 , J I = T* { I , J I «* Ct I,K)*TP2(K,J) 

30 CONTINUE 

I NTH ( 1 ,1 ) = NM 

INTH(2,1> a NVEC . „ 

CALL FECM<TP,EIG,VFC,HL,H,CNT,COLM,I NTH,NMAX) 

IF ( I NTH (1,1) «FQ. NM) GO 40 
PRINT 1000, lNTHUtl) 

STOP 

40 CONTINUE 

NORMALIZE VEC 
DO 50 J= 1 » NVE C 
SUMR = 0.0 
DO 45 I =1 1 NM 

45 SUMP = SUMR ♦ REAL C VEC( I, J ) )**2 * A I M AG ( V C C (I , J > » ** 2 
SUMR = SQRT(SUMR) 

OO 50 1=1, NM 

VEC (I ,J) = VEC ( I, J) /SUMR 
50 CONTINUE 

NORMALIZE AVFC 


DO 70 J=1 ,N VEC 
SUM = 0.0 

DO 60 1=1, NM % SUM 1 = 0.0 
DO 55 L =1 , NM 

55 SUM].= S UM1 +■ TP2(1,L»*VFC(L,JI 

60 SUM = SUM *■ AVEC( I, J >*SUM 1 
DO 7 0 I=1,NM 

AVFCt I , J t = A VEC ( I , J) /SUM 
70 CONTINUE 

IF ( IPRT .EC. 0) GO TO 100 
PRINT 2000, ( FIGt T ) ,1=1, NE IG ) 

IF ( IPRT .E Q. 1 ) GO TO 100 
PRINT 3000, ( ( VEC ( I , J I , J- 1 , NV FC 1 , 1 = 1, NM I 
PRINT 4000, { ( AVECt I , J) ,J = 1 ,NVEC) ,1 = 1, NM) 
100 CONTINUE 
RETURN 


1000 FORMAT <// , 3 1H 
2000 FORMAT l//,12H 
3000 FORMAT !//,13H 
4000 FORMAT!//, 24H 
C N0 


NUMBER OF FTGENVALU C S COMPUTED, 15) 
EIGENVALUES/ax,5(F14.3,E12.3m 
EIGENVECTORS,/, (IX ,5 ( F 14.3 , F 12 .3 ) ) ) 

ASSOC I AT FD El GEN VEC T ORS, / ,UX,5iE14.3, p 12.3))) 
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APPENDIX C 


USAGE DESCRIPTION OF LANGLEY LIBRARY 
SUBROUTINES USED BY PROGRAM MATCH 

Usage descriptions of the Langley library subroutines used by program MATCH 

are presented in this appendix. 

Langley Library Subroutine MATINV 

Language : FORTRAN 

Purpose: MATINV solves the matrix equation AX = B where A is a square coefficient 
matrix and B is a matrix of constant vectors. The solution to a set of simultaneous 
equations, the matrix inverse, and the determinant may be obtained. If the user does 
not want the inverse, use SIMEQ for savings in time and storage. For the deter- 
minant only, use DETEV, 

Use : CALL MATINV (A, N, B, M,DETERM, IPIVOT, INDEX, NMAX,ISCALE) 

A A two-dimensional array of the coefficients. On return to the calling 

program, A"* is stored in A. 

N The order of A; UNI NMAX. 

B A two-dimensional array of the constant vectors B. On return to calling 

program X is stored in B. 

M The number of column vectors in B. M = 0 signals that the subroutine 

is used solely for inversion, however, in the call statement an entry 
corresponding to B must still be present. 

DETERM Gives the value of the determinant by the following formula: 

DET(A) = (l0 100 ) ISCALE (DETERM) 
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I PIVOT A one -dimensional array of temporary storage used by the routine. 


INDEX A two-dimensional array of temporary storage used by the routine. 

NMAX The maximum order of A as stated in the dimension statement of the 

calling program. 

ISCALE A scale factor computed by the subroutine to keep the results of compu- 
tation within the floating point word size of the computer. 


Restrictions: Arrays A, B, IPIVOT, and INDEX are dimensioned with variable dimen- 
sions in the subroutine. The maximum size of these arrays must be specified in a 
DIMENSION statement of the calling program as: A (NMAX, NMAX), B (NMAX, M), 
IPIVOT (NMAX), INDEX (NMAX, 2). The orginal matrices, A and B, are destroyed. 
They must be saved by the user if there is further need for them. The determinant 
is set to zero for a singular matrix. 


Method: Jordan’s method is used to reduce a matrix A to the identity matrix I through a 
succession of elementary transformations: £ n , ^n-l» • • •» fi l* A = I. If these 

transformations are simultaneously applied to I and to a matrix B of constant vectors, 
the results are A -1 and X where AX = B. Each transformation is selected so that 
the largest element is used in the pivotal position. 

Accuracy : Total pivotal strategy is used to minimize the rounding errors; however, the 
accuracy of the final results depends upon how well-conditioned the original matrix is. 

Reference : Fox, L.: AN INTRODUCTION TO NUMERICAL. LINEAR ALGEBRA 

Storage : 542g locations. 
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Subroutine OPENMS 


Language : COMPASS 

Purpose : To open a random access file. 

Use : CALL OPENMS (U, IX, L,P) 
where 

U The logical unit number. 

IX The first word address of the index. 

L The length of the index. 

P P - 0 for numbered indexing. 

P = 1 for named indexing. 

Restrictions : OPENMS must be the first operation on a random access file. The file 
must be a disk file. For n index entries, the length of the index must be at least 
2n + 1 if using named indexing, whereas the index length must be at least n + 1 for 
numbered indexing. 

Method : OPENMS sets the first word in the index to a positive number for numbered 
indexing or to a negative number for named indexing. The random access bit, index 
address, and index length are set by OPENMS into the FET of the file for system 
communication. If the file already exists, the master index is read into central 
memory. 

Accuracy : Not applicable. 

References : None. 

Storage : 1030 locations. 
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Subprograms used : GETBA, SIO$, SYSTEM 

Error messages : (1) UNASSIGNED MEDIUM FILE XXXXXX 

(2) FILE DOES NOT RESIDE ON A RANDOM ACCESS DEVICE, 
XXXXXX 

(3) INDEX BUFFER IS OF INSUFFICIENT LENGTH, XXXXXX 


XXXXXX is the file name. Termination is abnormal in each case. 



APPENDIX C — Continued 
Subroutine WRITMS 

Language : COMPASS 

Purpose: To write a record on a random access file. 

Use : CALL WRITMS (U,FWA,N, I) 
where 

U The logical unit number. 

FWA The central memory address of the first word of the record. 

N The number of central memory words to be transferred. 

I The record number or record name depending upon the indexing mode 

set by the initial call to OPENMS. 

Restrictions : The file must have been opened by a call to OPENMS. 

Method : The specified record is written on the file and an address entered in the index 
to reference the record. 

Accuracy : Not applicable. 

References : None. 

Storage : 102s locations. 

Subprograms used : GETBA. SYSTEM, SIO$ 

Error messages : (1) UNASSIGNED MEDIUM, FILE XXXXXX 

(2) FILE WAS NOT OPENED BY A CALL TO SUBROUTINE OPENMS 

(3) INDEX BUFFER IS OF INSUFFICIENT LENGTH. 
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Subroutine RE AD MS 


Language : COMPASS 

Purpose : To read a record on a random access file. 

Use : CALL READMS (U,FWA,N,I) 
where 

U The logical unit number. 

FWA The central memory address of the first word of the record. 

N The number of words of the record to be transferred. 

I The record number or record name depending upon the indexing mode 

set by the initial call to OPENMS. 

Restrictions : The file must have been opened by a call to OPENMS. 

Method : The disk address of the record is determined by using the index. If n words 
are requested to be transferred and there are m words in the record, where m § n, 
m words are transferred. If m > n, n words are transferred. 

Accuracy : Not applicable. 


References : None. 

Storage : 13 lg locations. 

Subprograms used : GETBA, SYSTEM, SIO$ 
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Error messages: 


Termination is 


(1) UNASSIGNED MEDIUM, FILE XXXXXX 

(2) FILE WAS NOT OPENED BY A CALL TO SUBROUTINE OPENMS 

(3) RECORD NAME REFERRED TO IN CALL IS NOT IN THE FILE 
INDEX 

(4) *READ PARITY ERROR* 

(5) SPECIFIED INDEX IN THIS MASS STORAGE CALL .GT. MASTER 
INDEX OR IS ZERO, 

abnormal. 
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Subroutine EE CM 

Language : FORTRAN 

Purpose : To compute eigenvalues and eigenvectors of a complex N by N matrix. 

Use : CALL EECM (A, LAMBDA, VECT, HL, H, CNT, COLM, INTH, MAX) 

A A two-dimensional complex array of the input matrix. It is not 

destroyed. 

LAMBDA A one -dimensional complex array of eigenvalues. They are arranged 
in descending order of absolute magnitude. 

VECT A two-dimensional complex array of eigenvectors. Each vector is 

normalized so that the sum of the squares of the moduli of the 
components is unity. 

HL,H Two-dimensional complex temporary arrays. 

CNT, COLM One -dimensional complex temporary arrays. 

INTH A two-dimensional integer array. 

Upon entry - Before each CALL, set INTH as follows: 

INTH(1, 1) = N = order of matrix A. 

INTH(2, 1) = NV = number of eigenvectors to be computed. 

Upon return 

INTH(1, 1) = the actual number of eigenvalues computed. 

INTH(2, 1) is destroyed. 

MAX An integer, the maximum order of A. 

Restrictions : The calling program must type the following complex arrays and dimension 
them as follows: A(MAX, MAX), LAMBDA(MAX), VECT(MAX, NV), HL(MAX, MAX), 
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H(MAX, MAX), CNT(MAX), COLM(MAX). The integer array is dimensioned 
INTH(MAX, 2). 

Before each CALL to EECM, N and NV must be stored in the first 2 locations of 
INTH (see Use). 

The column dimension, NV, for VECT may be ^N. If no vectors are to be computed 
(INTH(2, 1) = 0), VECT need not be dimensioned, but it must appear as an argument 
in the call statement. 

The eigenvalues are not necessarily calculated in any absolute order, but are arranged 
in descending order of absolute magnitude prior to the calculation of the eigenvectors. 
Ten iterations per eigenvalue are allowed. In case of nonconvergence, the subroutine 
will return a value less than the order of the input matrix in INTH(1, 1). Thus, the 
user should test INTH(1, 1) upon return. If, then, it is less than the value of the 
number of vectors asked for, only that number of eigenvalues and eigenvectors is 
computed. If the number of eigenvalues computed is less than the order of the input 
matrix, the programer may want to use arbitrary shifts on the input matrix, or add 
a constant to the diagonal. Either change may eliminate the difficulty. Matrices apt 
to get nonconvergence are lower triangular with all equal eigenvalues, those with 
ones on the lower diagonal, and those with one as the Nth component of the first row 
and zeros elsewhere. 

If overflows or underflows occur, scaling the input matrix so that its largest element 
is in modulus about 1 will probably eliminate the difficulty. 

Equal computed eigenvalues return identical corresponding eigenvectors even though 
linearly independent vectors may exist. 

Method : The input matrix A is reduced to an upper Hessenberg matrix H by a sequence 
of elementary triangular and permutation matrices which make up a matrix P such 
that P“*AP = H. The QR algorithm is made use of in EECM by applying unitary 
similarity transformations to Hessenberg matrices, Hp Hj = P-^AP, H s = (hij(s)) 

- Qs T S) H s+1 - Q^H S Q S = Q^QsTgQs - TgQs where is the product of plane rota- 
tions, chosen so that T s is upper triangular. This process makes h^.j converge to 

zero and therefore h^ converges to an eigenvalue of A. When convergence is met 
( s') 

( h n, n-1 negligible), the Hessenberg matrix H s is deflated (i. e. , last row and column 
eliminated) and EECM proceeds with its leading principal submatrix (a new Hj) of 
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order one less. IE h^ s J l n _ 2 becomes negligible, the eigenvalues of the lower right- 
hand matrix of order two are calculated and EECM proceeds with the leading principal 
submatrix of order two less. It can be shown that convergence is accelerated by 
judiciously subtracting scalar matrices from the H s matrices. EECM actually 
replaces H s by H s - k s I so that k s is one of the eigenvalues p s or q s of the lower 
right-hand 2 x 2 matrix of H s . The choice of p s or q g is made on the basis of whether 
| h(s) _ p^j or |h^ - q g is a minimum. The shift technique is applied at each iteration. 

Two passes of the Wielandt inverse power method are used to calculate the eigen- 
vectors, Y i of H. Very little work is required for the second pass since the necessary 
elementary triangular and permutation matrices are stored in COLM and INTH(col. 2) 
(both internal storage areas). Finally, the eigenvectors of A, P Yj are calculated. 

The matrix P is in INTH (col. 1) and the lower part of H (internal arrays). 

The theory and a complete description of the algorithms appear in the first reference. 

Accuracy: The accuracy obtainable in computing the eigenvalues of input matrix A 
is usually related to the spectral radius, rho(A), of matrix A or more generally 
to some norm of A times the norm of its inverse. Hence, the greater 
rho(A) /min(abs(LAMBDA(l))), the fewer significant digits the smaller eigen- 
values may have. Accuracy also decreases as the order of the matrix 
increases. Close eigenvalues are usually less accurate than well separated ones. 

References: Wilkinson, J. H. : The Algebraic Eigenvalue Problem. Clarendon Press 
(Oxford), 1965. 

Householder, Alston Scott: The Theory of Matrices in Numerical Analysis. 
First ed., Blaisdell Pub. Co., 1964. 

Storage : 2745s locations. 

Subprograms used : None 

Timing: On Control Data 6000 computer, time for the actual solution of all eigenvalues 
and eigenvectors of a 30 by 30 matrix was 5. 2 seconds. This was about 5 times faster 
than routines presently in the Langley library. 
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